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CHAPTER  1:  INTRODUCTION 


RADAR  (Radio  Detection  And  Ranging)  is  an  electronic  device  for  measuring  the 
position  and  velocity  of  a  moving  object,  and  from  these  parameters  deduces  certain 
characteristics  of  that  object.  A  RADAR  operates  by  transmitting  an  electromagnetic  wave 
and  sensing  the  reflected  energy  in  space.  The  distance  or  range  to  the  object  from  the 
transmitter  is  determined  by  measuring  the  time  taken  by  the  pulse  to  travel  to  the  object 
and  back.  Since  the  electromagnetic  wave  travels  at  the  speed  of  light,  the  range  is  given 
by 


Range = 


(1-1) 


where  the  speed  of  light  c  =  3  TO8  m/sec  and  At  is  the  round  trip  travel  time  of  the  wave 
transmitted  and  reflected  back  to  the  source.  In  most  cases,  the  transmitted  wave  is 
periodic  and  the  time  elapsed  can  be  measured  from  the  origin  to  the  temporal  location  of 
the  peak  of  the  reflected  signal.  If  a  transmitted  pulse  is  received  after  the  second 
transmission,  the  measured  propagation  time  will  not  be  the  correct  one  since  the 
reference-transmitted  pulse  is  not  the  right  one.  This  situation  occurs  also  in  dwelling 
range  ambiguities.  The  maximum  range  without  ambiguity  would  depend  on  how 
frequently  the  transmission  occurs.  If  the  pulse  repetition  frequency  (PRF)  of  the 
transmitted  signal  is  low  and  the  time  interval  between  the  two  transmissions  is  long,  the 
maximum  range  that  can  be  measured  could  be  large.  However,  if  the  pulse  repetition 
frequency  (PRF)  of  the  transmitted  signal  is  high,  the  maximum  unambiguous  range  that 
can  be  measured  will  be  small.  Therefore,  the  maximum  range  rmax  without  any 
ambiguity  is  given  by 


bn  ax 


(1-2) 


where  fr  =  frequency  of  the  transmitted  signal.  If  the  target  is  located  beyond  this 
maximum  range,  the  system  will  predict  that  the  target  is  closer  than  the  actual  position 
due  to  folding  over  of  the  signal. 

The  velocity  of  the  target  can  be  detennined  from  the  change  in  the  carrier  frequency 
between  the  transmitted  and  received  signals  (Doppler  shift).  The  maximum  measurable 
velocity  without  any  ambiguity  also  depends  on  the  PRF  of  the  transmitted  pulse.  If  the 
Doppler  frequency  of  the  target  is  beyond  the  transmitted  PRF,  aliasing  may  occur  and  the 
real  velocity  of  the  target  cannot  be  obtained.  This  is  called  the  Doppler  ambiguity.  Here, 
the  term  “Doppler  frequency”  is  used  to  predict  the  velocity  of  the  target  from  the 
frequency  change  of  the  transmitted  pulse.  Obviously,  use  of  a  high  frequency  of  the 
transmitted  signal  could  work  for  faster  targets.  The  maximum  Doppler  frequency  without 
any  ambiguity  has  the  same  meaning  as  an  alias  free  sampling  used  in  the  Nyquist 
theorem.  That  is 


1 


(1-3) 


^max  /max^  0 


f,A{ 


where  Fmax  =  the  maximum  velocity  of  the  target, 

/max  =  maximum  Doppler  frequency, 

A  o  =  the  wavelength  corresponding  to  the  carrier  frequency  of  the  radar. 

In  other  words,  a  radar  operating  with  a  low  PRF  has  a  large  unambiguous  range  but  is 
ambiguous  in  Doppler  and  in  the  high  PRF  mode  there  is  ambiguity  in  range  but  not  in 
Doppler.  Therefore  the  trade-off  between  Doppler  ambiguity  and  range  ambiguity  is 
given  by 


(Maximum  range)(Maximum  Doppler  frequency)  = 


f  \ 

c 

(fA  o] 

cA  o 

U fr) 

l  2  2 

4 

(1-4) 


In  real  situations,  for  airborne  radars,  we  do  not  have  much  clutter  in  the  measurement 
of  Doppler  since  there  are  few  objects  moving  with  the  same  velocity  as  the  target. 
However,  more  clutter  would  be  found  in  the  measurement  of  range  and  that  uses  low 
PRF  rather  than  high  or  medium  PRF  radars.  To  improve  the  maximum  resolvable 
Doppler  frequency  (or  range),  many  techniques  have  been  proposed.  They  are  listed  as 
follows. 

•  Linear  carrier  FM  system  [  1 ,2] 

•  Sinusoidal  carrier  FM  system  [1,2] 

•  Use  of  a  Barker  coding  system  [2] 

•  Multiple  PRF  or  Staggered  PRF  systems  [  1 ,2] 

The  best  choice  depends  on  the  specific  application  and  on  the  choice  of  the  constraints. 
Generally  speaking,  a  multiple  PRF  system  performs  better  than  other  systems  [1].  Using 
several  fixed  PRFs  enables  one  to  discriminate  the  fold-overs  in  a  single  PRF  system  by 
comparing  the  responses  obtained  for  the  different  PRFs.  Thus  one  can  eliminate  either 
Doppler  or  range  ambiguities.  The  details  of  a  multiple  PRF  or  staggered  PRF  system  are 
described  in  the  later  sections. 

There  are  some  prior  methods  for  resolving  Doppler  and  range  ambiguities  in  a 
multiple  PRF  systems.  Ludloff  and  Minker  [3]  presented  a  curve  of  reliability  of  the 
velocity  measurement  from  simulations.  Vrana  [4]  dealt  with  the  problem  in  a  statistical 
manner.  The  two  steps  used  to  get  the  optimum  estimation  in  Vrana’ s  method  depends  on 
making  a  proper  decision  about  resolving  the  ambiguity  of  the  estimate  and  smoothing  of 
the  data. 

In  many  papers,  the  Chinese  remainder  (CR)  theorem,  which  will  be  explained  later  in 
the  report,  has  been  a  commonly  used  algorithm  to  resolve  Doppler  and  range  ambiguities 
[1,2],  In  spite  of  its  wide  usage,  the  CR  theorem  has  some  significant  defects  when 
applied  to  the  multiple  PRFs  systems.  First,  it  can  be  useful  when  there  is  a  single  target 
only.  If  there  are  several  targets  or  interferences,  which  have  the  same  received  power  at 
one  look  angle,  the  results  of  the  CR  theorem  would  be  ambiguous  since  there  will  be  too 
many  folded  Doppler  frequencies  to  be  resolved. 
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Moreover,  as  illustrated  in  Figure  1.1,  if  the  measured  frequency  fold-over  occurs  in 
the  region  of  the  main  lobe  clutter,  which  is  due  to  the  motion  of  the  aircraft  itself  with 
respect  to  the  ground,  the  resolution  of  the  target  is  not  clear.  Figure  1.1(a)  shows  a  target 
contact  and  the  clutter  due  to  the  ground.  The  various  parts  of  the  spectrum  would  be 
folded  due  to  the  PRFs  used  and  would  be  measured  as  in  Figures  1.1(b)  and  (c).  If  the 
image  falls  into  the  clutter  region,  as  shown  in  Figure  1-1  (c),  the  detection  cannot  be 
easily  performed.  Other  problems  with  this  approach  are  that  a  small  range  error  on  a 
single  PRF  can  cause  a  large  error  in  the  resolved  range  and  there  is  no  indication  that  this 
has  occurred. 


(a) 


Doppler  frequency 


Figure  1.1  Drawback  of  using  the  CR  theorem  (when  the  target  contact  folds  over  into  the 
clutter  region),  (a)  is  the  actual  contact  case,  (b)  folding  over  has  occurred  due  to  PRF1  in 
which  case  the  measurement  is  not  in  the  clutter  region  and  (c)  folding  over  has  occurred 
due  to  PRF2  where  the  measurement  is  in  the  clutter  region. 
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To  overcome  this  problem  due  to  the  use  of  the  CR  theorem,  a  clustering  algorithm  has 
been  proposed  by  Trunk  and  Rockett  [5].  They  used  a  variance  for  each  PRF  data  to  test 
if  it  is  the  real  Doppler  or  range.  This  method  will  be  explained  in  detail  in  section  3.2. 
Even  with  the  clustering  algorithm,  the  fundamental  drawback  of  the  CR  theorem  based 
approach  cannot  be  resolved  since  the  solution  is  obtained  through  numerical  techniques. 
However,  none  of  the  methods  deals  with  the  problem  of  obtaining  a  single  estimate  from 
the  plethora  of  unevenly  spaced  data  obtained  for  a  multiple  PRF  RADAR  system.  As 
illustrated  in  Figure  1.2,  the  previous  researches  separate  the  solution  space  for  each  PRF 
and  consider  it  as  a  combination  of  single  PRF  systems.  In  the  frequency  domain,  Figure 
1.2(d,  e,  f)  has  a  common  peak  at  around  lOOrad/sec.  Instead  of  dealing  with  one  PRF  at  a 
time,  one  can  think  of  assembling  the  data  from  all  the  PRFs  simultaneously  in  the  time 
domain  as  shown  in  Figure  1 .3(a).  The  problem  of  Figure  1 .3(a)  is  then  reduced  to  finding 
the  spectrum  of  an  unevenly  spaced  sampled  signal.  If  we  can  find  the  spectrum  of  a 
nonuniformly  sampled  signal,  we  would  not  have  the  problem  of  having  to  solve  for  the 
congruence  and  deal  with  the  various  disadvantages  due  to  the  CR  theorem.  Moreover, 
the  system  could  deal  with  multiple  targets  and  it  could  also  resolve  the  blind  speed  and 
the  blind  range  problems  with  the  removal  of  ambiguities. 


(a)  Signal  sampled  by  PRF  1 


(d)  Frequency  response  of  (a) 


Figure  1.2  Sampled  signal  at  multiple  PRFs  and  their  frequency  domain  response  (DFT).  Signal  frequency 
was  100  rad/sec  and  all  the  frequency  domain  results  are  aliased.  Observe  that  (d,  e,  f)  has  a  common  peak 
around  100  rad/sec  which  is  the  real  frequency. 
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(a)  Signal  sampled  by  PRF  1,  2  and  3  simultaneously 

- > - t - . -  1  o 


(b)  Frequency  response  of  (a) 


Figure  1.3  Unevenly  spaced  signal  and  its  spectrum  (DFT). 

In  this  research,  we  want  to  obtain  a  frequency  domain  response  from  a  nonuniformly 
sampled  sequence.  Seven  methods  for  unevenly  spaced  data  analysis  have  been  studied  in 
this  report  and  summarized.  They  have  also  been  used  to  simulate  results  for  a  multiple 
PRF  case.  The  seven  algorithms  consist  of 

•  Polynomial  interpolation  (Lagrange  and  Cauchy  type) 

•  Chinese  remainder  theorem  and  clustering  algorithm 

•  Least  squares  curve  fitting  of  a  complex  sequence 

•  Multi-resolution  (Quadrature  mirror  filter)  analysis 

•  Iterative  method 

•  Orthogonal  expansion  by  a  set  of  polynomials  (Legendre  polynomials  and 
Hermite  polynomials) 

•  Estimation  of  an  analog  frequency 

Some  of  these  methods  generate  an  evenly  spaced  sequence  from  the  unevenly  spaced 
data.  Flence,  for  those  methods,  the  FFT  is  then  utilized  to  estimate  the  frequency 
components  from  the  evenly  spaced  sequences.  The  matrix  pencil  method  [7,  8]  can  also 
be  used  to  efficiently  extract  the  parameters  with  a  higher  resolution  of  the  frequency 
domain  sequence  obtained  from  the  FFT.  The  Matrix  Pencil  Method  has  been  discussed 
in  Appendix  B.  In  Chapter  3,  all  the  above  approaches  have  been  described.  Computer 
simulated  examples  have  been  presented  for  all  of  them.  The  results  have  been  compared 
to  investigate  which  approach  is  applicable  to  the  radar  application.  A  summary  of  the 
methods  is  given  in  Chapter  4,  which  is  also  the  conclusion. 

Additional  benefit  of  using  a  direct  analysis  of  unevenly  spaced  data  is  that  it  can 
reduce  the  distortion  in  the  spectrum  of  a  signal  affected  by  noise  due  to  the  correlation 
associated  with  each  of  the  frequency  component  [11].  It  is  generally  known  that  if  the 
sampling  is  completely  random,  and  is  a  Poisson  process  [9],  then  the  spectrum  of  that 
sequence  is  alias  free.  A  proof  of  this  statement  is  given  in  the  references  [12,  13]  and  a 
sketch  of  it  is  given  in  Appendix  A. 
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CHAPTER  2:  DESCRIPTION  OF  MULTIPLE  PRFS  SYSTEM 

In  this  section,  the  performances  of  multiple  PRF  systems  are  investigated.  It  will  be 
seen  that  they  have  an  enhanced  performance  when  compared  to  that  of  a  single  PRF 
system.  First,  consider  an  ambiguity  function  (AF)  which  is  a  tool  often  used  for 
characterizing  ambiguities.  We  assumed  a  transmitted  radar  signal  of  the  general  form 

g0(f)  =a(f )cos  [27rFct  +  x¥(t)\,  (2-1) 

where  a(t)  is  the  envelope  of  the  signal;  Fc  is  the  carrier  frequency  and  vF(t)  is  the 
phase.  If  this  signal  illuminates  a  target  moving  at  a  speed  v,  the  transmitted  signal 
undergoes  a  frequency  shift  due  to  the  Doppler  effect  and  the  mathematical  expression  for 
the  received  echo  becomes 

g(t,a)=a(ta)  cos  [2^Ffta  +  vP(t«)],  (2-2) 

where  a  is  a  scale  factor  controlled  by  the  Doppler  effect  and  is  given  by 

c  -  v 

a  = - ,  (2-3) 

c  +  v 

where  c  =  velocity  of  light.  The  Doppler  frequency  of  the  target  is 

Fd={l-a)Fc,  (2-4) 

If  we  pass  gif,  a)  through  a  signal  conditioner,  which  converts  the  carrier  frequency  to  the 
intennediate  frequency  (IF)  radar  signal  and  then  through  a  low  pass  filter,  L,  to  get  the 
base  band  signal,  one  obtains 

g(t,a)  =  L[^r-'g(t,a)]^c(al)  '»  (2-5) 

Assume  that  the  phase  function  vF(t)  is  zero  and  the  envelope  is  a  flat  topped  pulse,  so 

T  / 

that  a  (a  t)  =  1  and  for  |r|  <  py  where  Tp  is  width  of  the  pulse  in  a  period,  (2-5)  will 
become 

g(t,a)  =  e2Mi('~a)Fc‘ ;  |,  -  kT\  <  Ty2  ,  (2-6) 

where  T  is  the  period  of  the  base-band  pulse  and  k  is  an  integer. 

The  complex  ambiguity  function  (CAF)  of  g(t)  is  defined  through 
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(2-7) 


A (t,a)=  J  g(u-t, l)  g{u, a)  du. 

—oo 

A 

If  a  =1  in  g(u  -  /  ,1)  this  implies  that  there  is  no  Doppler  shift  and  the  waveform  is  shifted 
by  t  along  the  u  axis.  Next,  this  is  multiplied  by  the  original  signal  and  integrated  with 
respect  to  u.  This  results  in  the  familiar  form  of  an  auto-correlation  function  for  each 
Doppler  frequency.  If  there  is  no  Doppler  shift,  then  the  CAF  becomes  merely  an  auto¬ 
correlation  function  and  can  be  used  to  measure  the  unambiguous  range.  Actually,  the 
main  lobe  of  the  CAF  measures  the  range  ambiguity  of  a  signal  when  there  is  no  Doppler 
shift.  When  t  =  0  this  implies  that  there  are  no  range  shifts  and  the  unambiguous  Doppler 
can  be  determined  from  the  main  lobe  of  the  CAF.  Therefore  the  CAF  measures  the 
maximum  of  the  Doppler  shift  and  range  which  can  be  resolved  by  a  given  signal  model. 
The  phrase  “unambiguous  function”  is  more  suitable  instead  of  the  ambiguity  function 
since  a  large  value  of  CAF  implies  a  larger  domain  for  the  Doppler  and  the  range 
estimation  without  ambiguity. 


Figure  2.1.  A  baseband  rectangular  pulse  train. 

Figure  2. 1  is  a  baseband  time  domain  signal  due  to  a  rectangular  pulse  train  a{t) .  The 
width  of  the  pulse  T  is  0.2msec  and  the  period  T  is  0.5msec  (2kHz).  The  transmitted 
signal  is  modulated  by  cos  {2nFct) .  The  numerical  integration  of  (2-7)  using  (2-4)  and 
(2-6)  is  shown  in  the  Figure  2.2(a).  As  a  result, 


A(c«)=Z 

k 


(  T 
kT+  ”/ 


J  e2mF“udu 


\kT~  72 


(2-8) 
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Figure  2.2(b)  provides  the  corresponding  contour  plot.  The  width  of  the  main  lobe  is 
determined  by  the  duration  of  the  baseband  pulse. 


(a) 


x  10‘3 


Figure  2.2  (a)  3  dimensional  plot  of  CAF  for  the  single  PRF  case  (b)  Contour  plot  of 

CAF. 

For  a  multiple  PRF’s  CAF,  only  the  envelope  term  will  change  and  the  summation  in 
(2-8)  shall  be  computed  for  different  T p  s.  Figure  2.3  is  the  baseband  envelope  for  a  two 
PRFs  system  where  1kHz  and  1.5kHz  have  been  chosen  as  the  two  PRFs  which  have  the 


same  number  of  pulses  in  time  as  in  Figure  2.1  for  comparison  with  a  single  PRF  case.  (2- 
7)  will  become 


A(^«)  =  Z 


X  fe2mF"udu 


\ll  1 


(2-9) 


Figure  2.3  A  multiple  PRF  (1kHz  and  1.5kHz)  baseband  pulse. 

Note  that  the  pulse  exists  for  tn  <  t  <  tl2  and  /  goes  up  to  4  since  it  is  repeated  after  every 

4  pulses.  Figures  2.4(a)  and  (b)  are  the  corresponding  plots  of  the  contour.  By  comparing 
Figure  2.2  with  Figure  2.4,  it  is  seen  that  a  2  PRF  system  performs  better  than  the  single 
PRF  system  since  it  has  a  larger  unambiguous  region  and  the  width  of  the  main  lobe  is 
wider  than  that  of  the  single  PRF  case.  One  can  observe  from  Figure  2.1  and  2.4  that  the 
maximum  unambiguous  range  for  the  2  PRFs  case  has  increased  to  2.0msec  which  is  4 
times  that  of  the  single  PRF  case  of  0.5msec.  Therefore  a  multiple  PRF  system  would 
have  wider  unambiguous  regions. 

The  maximum  resolvable  frequency  and  range  in  a  multiple  PRF  system  is  determined 
from  the  set  of  PRFs.  Consider  the  3  PRFs  ( PRFl,PRF2  and  PRF, )  which  sample  the 

radar  signal  and  compare  its  performance  to  the  single  PRF  ( PRF2 )  system.  First  assume 
that  the  3  PRFs  are  relatively  prime  numbers  with  respect  to  each  other.  The  maximum 
unambiguous  Doppler  increases  with  the  product  of  the  PRFs  since  the  frequency 
spectrum  from  each  PRF  will  coalesce  at  a  frequency  multiple  of  those  PRFs.  That  is, 
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ID  = 


(2-10) 


PR  I]  ■  PRF2  ■  PRF 3 

prf2 

where  ID  is  the  maximum  Doppler  increment. 


i 
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Figure  2.4  (a)  3  dimensional  plot  of  CAF  for  a  2  PRF  system  and  (b)  Contour  plot  of  CAF 
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(b) 
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As  an  example,  the  maximum  resolvable  Doppler  frequency  of  a  4Hz  PRF  system  is 
%  Hz  from  (1-3)  and  the  maximum  resolvable  Doppler  frequency  of  a  3  PRF  system 

consisting  of  3,  4  and  5Hz  is  3  -4  -5/^  Hz.  If  the  PRFs  are  not  relatively  prime  numbers, 

the  maximum  discemable  Doppler  would  be  the  Least  common  multiplier  ( l.c.m )  of  those 
PRFs.  That  is,  it  will  be  increased  by  the  factor  ID,  where 

ID  =  IjCjnjPRFy ,  PRF2  ,...,PRFn  ) 

PRF2 


Hence,  the  performance  would  be  ID  times  better  compared  to  that  of  a  single  PRF 
system,  which  is  PRF2 .  Typically  the  performance  is  enhanced  by  the  maximum  Doppler 
increment. 

The  maximum  resolvable  range  is  also  increased.  In  terms  of  PRF,  i.e.,  the  range  is 
proportional  to  Yprf  anc*  range  increment  is  proportional  to  the  PRF.  When  the 


PRFs  are  relatively  prime,  the  period  of  the  repeating  pulse  will  be  the  l.c.m  of 


1/ 

/  PRFl  ’ 


1/ 

/PRF. 


and  \ 


PRF \ 


which  is  unity.  The  perfonnance  is  enhanced  by 


IR  = 


7 - =  PRF,. 

prf. 


(2-12) 


where  IR  is  the  maximum  range  increment. 

It  is  important  to  note  that  it  does  not  matter  how  many  PRFs  exist  but  the 
enhancement  in  range  resolution  given  by  (2-12)  occurs  when  the  PRFs  are  relatively 
prime  because 


l.c.m\ 


PRF,  ’  7 PRF. ’  7  PRF, 


=  l.c.ml 


/ PRF ;  ’  7  PRF, 


=  1. 


(2-13) 


As  an  example,  consider  the  maximum  resolvable  range  of  the  4Hz  system  which  in 
meters  is  4)’  obtained  from  (1-2)  and  the  maximum  resolvable  range  of  the  3,  4  and 

5Hz  PRF  system  are  meters.  If  the  PRFs  are  not  relatively  prime  numbers,  then 


IR  = 


prf2 

g.c.d{PRFx ,  PRF, ,  PRF/  ’ 


(2-14) 


where  g.c.d  =  greatest  common  divisor.  Then  the  maximum  resolvable  range  is  the 
maximum  range  increment  times  that  of  a  single  PRF  system  of  value  PRF, .  The  product 
of  (2-11)  and  (2-14)  results  in  a  performance  enhancement  of 
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l.c.m  (PRF, ,  PRFj ,  PRF, 

ID  IR= - 7 - - - - - ^4 

g.c.d  {PRF] ,  PRF, ,  PRF, 


(2-15) 


times  that  from  a  single  PRF  radar  system. 


CHAPTER  3:  DESCRIPTION  OF  THE  VARIOUS  METHODS 

In  this  chapter,  various  methods  to  obtain  the  spectrum  from  a  set  of  nonuniformly 
sampled  data  are  described.  The  first  class  of  methods  presented  process  the  nonunifonnly 
spaced  data  through  a  spatial  (or  time)  domain  interpolation  to  a  uniformly  sampled  case 
and  then  uses  the  conventional  FFT  or  DFT  to  get  the  frequency  response.  The  second 
class  of  methods  directly  obtains  the  spectrum  from  a  set  of  nonunifonnly  sampled  data. 
The  basic  polynomial  interpolation  and  the  iterative  methods  fall  in  the  first  category,  and 
the  rest  of  the  methods  described  in  this  chapter  belong  to  the  second  group.  Usually,  the 
approach  using  spatial  domain  interpolation  requires  much  more  densely  spaced  data 
samples  which  have  a  greater  sampling  rate  than  that  of  the  Nyquist  sampling  rate  to 
provide  meaningful  results  for  the  spectrum. 

Note  that  the  Nyquist  sampling  rate  for  a  nonuniformly  sampled  data  can  be  defined  as 
the  average  sampling  rate  of  the  sequence  and  it  can  be  shown  that  if  the  average  sampling 
rate  exceeds  twice  the  maximum  frequency  of  the  actual  signal,  then  the  signal  can  be 
perfectly  recovered  [12,  13].  This  is  an  extension  of  the  actual  Nyquist  sampling  theorem 
corresponding  to  the  nonunifonnly  sampled  data  case. 


3.1  INTERPOLATING  IN  SPATIAL  DOMAIN  BY  POLYNOMIALS 

3.1.1  Lagrange  Interpolation  Polynomials 

There  are  some  obvious  ways  to  generate  uniformly  spaced  data  from  a  nonuniformly 
sampled  sequence.  Interpolation  using  polynomials  is  one  way  to  do  it.  One  of  the 
simplest  and  direct  interpolation  schemes  involves  the  use  of  the  Lagrange  interpolation 
polynomial  which  fits  a  set  of  N  data  points  by  an  (N- 1  )th  degree  polynomial  where  N  is 
the  given  number  of  data  samples.  The  interpolation  may  become  smooth  when  there  are 
enough  data  points  in  one  period  of  the  signal.  The  results  of  which  would  be  quite 
acceptable  if  the  spacing  is  not  very  random  (small  deviation  form  uniform  spacing). 
However,  this  approximation  merely  provides  a  base  line  of  the  interpolation  and  does  not 
exploit  any  property  from  the  frequency  domain  like  the  signal  is  periodic  in  nature. 
Previous  research  has  indicated  that  if  the  reconstruction  of  the  signal  from  random 
samples  is  performed  through  the  use  of  interpolation  by  polynomials,  the  errors  in  the 
reconstruction  are  acceptable  only  if  the  sampled  frequency  of  the  signal  components  is  at 
least  four  times  than  that  of  the  Nyquist  sampling  rate  [14,  15].  Once  we  get  the  uniformly 
spaced  data,  the  Fast  Fourier  Transformation  (FFT)  or  the  Matrix  pencil  method  can  be 
used  to  estimate  the  frequency  domain  parameters. 

The  Lagrange  interpolation  formula  between  the  sample  points  evaluates  the  function 
through  the  following  interpolation 
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,(3-1-1) 


n,  ,  (x-x2)  (x-x3)  —  (x-xN) 

P(x)  =  7 - w - - 7/(Xl }  + 

(*i  -x2)  (x1  -x3)---(x1  -x*) 

(x  —  Xj )  (x-x3)---(x-xJV) 


(x2  -Xj)  (x2  -x3)---(x2  -Xw) 


/(x2  )  -t - + 


(x  —  Xj  )  (x  —  x2  )  (x  —  x 2V_1 ) 
(Xjv-Xj)  (xN-x2)  —  (xN-xff_l) 


f(xN) 


where  xk  is  the  location  of  the  nonuniformly  spaced  sampled  points  k  =  1,  N,  and. 

/(xA.)  is  the  value  of  the  signal  at  xk .  This  amounts  to  fitting  a  (/V- 1  )-th  degree  polynomial 

through  the  N  data  points.  One  practical  problem  associated  with  this  technique  is  that  when 
the  number  of  data  points  becomes  large,  Equation  (3-1-1)  has  to  deal  with  very  large 
numbers  because  there  are  coefficients  with  A-th  power  of  the  time  argument.  This  can  cause 
numerically  unstable  results.  To  prevent  the  order  of  the  polynomials  from  being  a  very  large 
number,  a  time  domain  scaling  can  be  performed.  The  value  of  the  function  is  scaled  between 
-1  to  1,  while  also  checking  for  N  so  that  it  does  not  become  a  large  number.  More 
consideration  should  be  given  to  the  edges  of  the  data  sequence  since  at  those  regions  the 
polynomial  may  not  accurately  fit  the  data. 

As  an  example  consider  the  following  signal 

f(xk )  =  e2'lraki  +  2e~3,5'2*v  -  2.5e4'2mk‘  (3-1-2) 

andx  is  a  randomly  generated  number  between  -1  to  1.  41  samples  have  been  chosen  to  make 
the  average  sampling  time  to  be  0.05  (average  sampling  frequency  =  20Hz).  Note  that  the 
maximum  frequency  of  the  signal  is  4Hz.  It  is  one-fifth  that  of  the  Nyquist  sampling  rate  since 
the  signal  is  complex  and  in  which  case  the  minimum  sampling  rate  for  the  perfect  recovery 
of  the  signal  is  equal  to  the  maximum  frequency  of  the  signal  so  as  to  completely  eliminate 
the  ambiguities  in  the  results. 

Note  that  the  spectrum  of  ejcot  exists  only  along  the  positive  axis  if  co  >  0 ,  while  the 
spectrum  of  sinfin  /)  exists  for  both  positive  and  negative  frequencies.  The  various  time 

domain  signals  are  shown  in  Figure  3.1.1  along  with  the  original  signal.  The  corresponding 
spectrum  is  given  in  Figure  3.1.2  along  with  the  results  for  the  unifonnly  sampled  data  using 
FFT.  Decreasing  the  sampling  rate  or  increasing  the  signal  frequency  causes  the  interpolation 
to  become  inaccurate.  There  is  no  guarantee  of  convergence  of  this  process  to  the  original 
signal  unless  infinite  samples  in  one  period  are  taken.  Generally,  these  interpolations  perform 
poorly  in  the  computation  of  the  spectrum  for  a  nonuniformly  spaced  signal  as  compared  to 
other  methods  described  in  this  report.  However,  this  method  offers  a  base  line  comparison  in 
the  analysis  of  nonuniformly  spaced  data  when  the  sampling  rate  is  much  less  than  the 
Nyquist  rate.  In  addition,  it  can  also  be  easily  implemented  in  hardware. 

3.1.2  Cauchy’s  Method 

Cauchy’s  method  is  a  technique  of  finding  a  rational  polynomial  which  will  fit  a  given  data 
sequence.  Previous  researches  have  successfully  interpolated  data  from  an  electromagnetic 
system  using  this  approach  [16,  17].  A  brief  introduction  and  derivation  of  the  Cauchy’s 
method  and  an  application  to  the  nonuniformly  sampled  interpolation  are  presented  in  this 
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section.  The  main  objective  of  the  Cauchy’s  method  is  to  find  the  coefficients  and  the  orders 
of  the  polynomials  for  the  numerator  and  the  denominator. 

Assume  that  the  signal  can  be  approximated  by  the  rational  polynomial 


P 


k= 0 


(3-1-3) 
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Figure  3.1.1  Use  of  the  Lagrange  interpolation  polynomial  (magnitude  only)  to  a  time  domain  data. 


Frequency  [Rad/sec] 

Figure  3.1.2  FFT  of  the  time  domain  interpolated  data  of  Figure  3.1.1  due  to  Lagrange  interpolation 
(magnitude  only) 
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Then  the  Cauchy  problem  can  be  stated  as:  Given  //(x,)  for  i=l,...,N,  find  p,  q, 
ak(k=0,...,p)  and  bk(k=0,...,q). 

By  enforcing  the  equality  of  both  sides  in  (3-1-3),  the  result  is  obtained  as 

A{xi)=H(xi)-B(xi)  (3-1-4) 

or  equivalently 

a0  +  a]xj  H - 1 -apxf  —  i/(x,.)  b0  -  H{xt)  blxi - //(x,)  bqxqi  =0  for  i=l,  ...,7V.  (3-1-5) 

A  matrix  fonn  of  this  equation  will  then  become 

Amata  =  Bma,b  (3-1-6) 


1  Xi 

•• 

Cl\ 

h 

A 

1  X2 

••  x2 

D 

?  °mat  ~ 

H(x  2)  H{x  2)^2 

■  H(xl)x2 

a2 

h  - 

b2 

^ mat 

J  xN  • 

'■  XN_ 

H(xn)  H(xn)xn  ■ 

■  h(xn)xn_ 

,  Cl  — 

ap 

?  °  ~ 

or 

“  =»  (3-1-7) 

b 

The  singular  values  of  [Amat  \  -Bmat  ]  can  be  obtained  by  using  the  singular  value 
decomposition.  The  number  of  nonzero  singular  values  of  [ Amat  \  -Bmat\  will  be  the  sum 

of  order  of  the  denominator  and  the  numerator.  It  provides  some  guidance  in  estimating 
the  values  of  p  and  q.  If  z  is  the  number  of  nonzero  singular  values,  then  p  and  q  should 
satisfy  the  relationship 

z  =  p  +  q  +  2  (3-1-8) 

and  p  is  chosen  such  that  q=p+l.  To  obtain  a  and  b,  apply  a  QR  decomposition  to  Amat . 

That  is, 

QRa-Bb  =  0  (3-1-9) 

Since  Q  is  an  orthogonal  matrix  (Q  1  =  Q  '  and  therefore 

Ra-QTBb  =  0  (3-1-10) 
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The  rank  of  R  is  determined  by  the  order  of  the  numerator  polynomial  which  is  p+1 ,  and 
(3-1-10)  becomes 


R-QTB 


o 

a 

1 

<N 

i _ 

a 

u 

b 

L°  r2li 

b 

0 

(3-1-11) 


where 


value  decomposition  of  r22 ,  i.e., 


Therefore  b  can  be  obtained  from  the  singular 


r22b  =  UYVTb  =  0  (3-1-12) 

From  the  theory  of  total  least  square  (TLS)  [18],  the  solution  to  (3-1-12)  is  the  last  column 
of  the  matrix  V. 


b  =  (3-1-13) 

Therefore 

a  =  RlQTBb  (3-1-14) 

The  same  signal  presented  in  the  previous  section  has  been  used  as  an  example.  That  is, 

f(xk )  =  e2'lnXki  +  2e~3'5'2wCti  -  2.5e4'2mki 

and  x  is  a  randomly  generated  number  between  -1  to  1.  41  samples  of  the  data  have  been 
chosen  to  achieve  the  average  sampling  time  of  0.05  (average  sampling  frequency=20Hz). 
The  time  domain  result  is  shown  in  Figure  3.1.3  and  has  been  compared  to  that  of  the 
original  signal.  The  corresponding  spectrum  of  the  signal  is  given  in  Figure  3.1.4  along 
with  the  FFT  of  an  evenly  spaced  sequence.  The  orders  of  the  polynomials  are  chosen  to 
be  q=31  and  p=30. 
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Original  Signal 

. 

Sampled  Signal 

o 

Uniformly  Interpolated  Signal 

Figure  3.1.3  Time  domain  interpolation  using  the  Cauchy’s  method,  (magnitude  only). 


Frequency  [Rad/sec] 

Figure  3.1.4  FFT  of  the  waveform  interpolated  by  the  Cauchy’s  method  as  shown  in  Figure  3.1.3. 
(magnitude  only). 


3.2  CHINESE  REMAINDER  THEOREM  AND  THE  CLUSTERING  ALGORITHM 

The  Chinese  remainder  theorem  and  the  clustering  algorithm  is  a  different  approach  as 
compared  to  the  other  methods  described  in  this  chapter.  In  this  case,  one  estimates  the 
target  Doppler  frequency  from  a  set  of  frequencies  computed  from  each  PRF  by  taking  the 
conventional  FFT  of  the  evenly  sampled  data.  Since  we  want  the  maximum  Doppler 
frequency  of  the  signal,  which  exceeds  the  Nyquist  sampling  rate,  aliasing  may  occur.  If 


17 


the  desired  frequency  component  is  larger  than  the  sampling  frequency,  then  a  fold  over 
along  the  sampling  frequency  will  occur. 

Therefore,  the  aliased  frequency  is  the  residue  of  the  result  of  the  division  of  the 
original  signal  frequency  by  various  integers.  This  is  shown  in  Figure  3.2.1.  Here,  the 
original  frequency  of  the  signal  is  65kHz  and  it  is  aliased  when  sampled  at  a  rate  of  15kHz 
and  is  measured  at  5kHz.  Figure  3.2.1(b)  illustrates  the  results  for  the  same  problem  when 
the  sampling  frequency  is  20kHz  and  the  spectrum  is  still  measured  at  5kHz.  Figure 
3.2.1-(c)  shows  the  case  when  the  sampling  frequency  is  25kHz  with  the  measured  aliased 
value  at  15kHz.  This  problem  can  be  solved  using  the  following  equations 

65 kHz  =  m ,  x  1 5 kHz  +  5kHz  , 

65 kHz  =  m2  x  20 kHz  +  5 kHz ,  (3-2-1) 

65 kHz  =  m3  x  25 kHz  +  15 kHz  . 

One  can  easily  see  that  mx=  4,  m,=3  and  m,=2.  If  the  information  exists  only  for  the 
sampling  frequencies  and  the  measured  aliased  frequencies  and  the  signal  frequency  is 
unknown,  the  calculation  for  the  minimum  value  of  the  signal  frequency  is  done  by 
choosing  the  minimum  multiple  of  the  integers  ml ,  m1  and  m3  ■ 


(a)  Sampling  frequency  =  15kHz 


Meas 

jred  frequency  1 

- 

0  5  10  15  20  25  30 

frequency  [kHz] 

(b)  Sampling  frequency  =  20kHz 
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- 
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(c)  Sampling  frequency  =  25kHz 
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d  frequency  3 
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Figure  3.2.1  Results  for  Doppler  ambiguity 
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The  equation  still  holds  for  any  integer  multiples  of  m] ,  m2  and  m3-  In  this  case,  the 
maximum  resolvable  Doppler  frequency  would  be  the  Least  common  multiplier  ( l.c.m )  of 
those  sampling  frequencies.  This  is  computed  from, 

fo  =  PRfV>h  +/,,/„  =  PRF2m2  +f2,fa  =  PRF3m3  +  f3 .  (3-2-2) 

where  /0  =  target  Doppler  and  ml,m2  and  m3  are  integers,  or 

fo  =  fi  mod (PRFl ) ,  f0=f2  mod (PRF2 ) ,  f0  =  f3  mod (PRF3 ) ,  (3-2-3) 

where  mod  is  defined  through  the  expression  a  mod(h)  =  nb  +  a  and  n  =  an  integer.  The 
problem  then  becomes  one  of  solving  a  set  of  modulus  equations,  so  called  simultaneous 
congruences,  which  should  be  satisfied  at  the  same  time. 
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Figure  3.2.2  Multiple  PRF  resolves  range  ambiguity. 


The  same  procedure  can  be  applied  to  determine  the  range  as  shown  in  Figure  3.2.2. 
The  measured  distances  for  each  PRF  are  denoted  by  x1 ,  x2  and  ,v3 .  The  real  range  would 
be 


x0  =  Tlnl  +  Xj ,  xD  =  T2n2  +  x2 ,  xa  =  T2n2  +  x 


(3-2-4) 
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where  x0  =  target  range,  and  nl,n2  and  n3  are  integers,  with  7]  =  YpRI 71  ’  ^  =  YpRF2  ’ 
T3  =  YpRF3  ■  Equivalently 


x0  =  x,  mod(lj),  x0  =  x2  mod(r2),  x0  =  x3  mod(r3) .  (3-2-5) 

Note  that  one  can  measure  the  time  instead  of  the  distance  since  the  distance  to  the  target 
equals  c  . 

The  most  common  algorithm  to  resolve  a  set  of  simultaneous  congruences  like  in  this 
case  is  the  Chinese  remainder  theorem.  This  is  the  most  common  technique  used  currently 
and  much  research  has  been  done  establishing  the  credibility  of  this  approach  [3-5].  First, 
consider  a  case  of  two  congruences 


From  (3-2 -6-a), 


x  =  b  mod(n), 

(3-2-6-a) 

x  =  a  mod(m). 

(3-2-6-b) 

x  =  b  +  nt 

(3-2-7) 

and  from  the  second  equation  one  observes  that  t  must  satisfy  the  condition 

a  +  mt  =  b  modi/?)  (3-2-8) 

mt  =  (b  -  o)mod(n) .  (3-2-9) 


or 


According  to  the  general  rules  just  derived,  the  linear  congruence  in  t  can  only  have  a 
solution  when  the  greatest  common  divisor  denoted  by  g.c.d(m,  n )  can  be  divided  by  b-a. 
When  this  is  the  case  the  congruence  (3-2-9)  may  be  divided  by  d  and 


m  b-a  , 

—  t  = - mod 

d  d 


f  n  3 
ydj 


(3-2-10) 


Let  t0  be  some  particular  solution  of  this  congruence  and 


x0  =  a  +  mt0 , 


which  is  a  solution  of  (3-2-9).  The  general  solution  of  (3-2-9)  is  then 


t  =  t0  mod 


f  n  3 
\dj 


(3-2-11) 


so  that  it  can  be  written  as 
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(3-2-12) 


n 

t  =  t0+u~, 

where  u  is  some  integer.  The  resulting  general  solution  of  the  original  congruence  is 

^  n  3  mn 

x  =  a  +  m  ta  — u  =  xn  +u - 

d  )  0  d 

or 

x  -  x0  mod  [l.c.m(m,  n)\ , 

since  l.c.m  ( m,n )  =  is  the  Least  common  multiplier  for  m  and  n. 

When  one  considers  a  set  of  algebraic  congruences  for  several  modulus  and  x0  is  the 
number  satisfying  all  of  them,  it  is  clear  that  if  one  adds  any  multiples  of  l.c.m  of  all 
modulus  of  x0 ,  the  resulting  number  will  also  be  a  solution.  Therefore,  with  several 
possible  moduli  it  is  apparent  that  the  number  of  different  solutions  is  given  by  the 
incongruent  solution  corresponding  to  the  l.c.m  of  the  various  modulii. 

When  several  simultaneous  congruences  are  given 

x  =  a]  mod(/«,) ,  x  =  a2  mod(m2) ,  x  =  a3  mod(/«3) ,  (3-2-15) 

then  the  solution  can  be  found  by  repeated  application  of  the  method  given  above.  One 
combines  the  first  two  congruences  and  finds  a  single  congruence  as 

x  =  a0  mod  [l.c.m(mx ,  m2  )],  (3-2-16) 

which  can  replace  (3-2-10).  This  in  turn  is  solved  in  conjunction  with  the  third,  and  so  on. 
One  sees  that  if  there  exists  a  solution  of  the  congruences  (3-2-10),  then  there  is  only  a 
single  solution,  with  respect  to  a  modulus  that  is  the  l.c.m  of  all  the  modulus  mi . 

The  necessary  and  sufficient  condition  for  a  set  of  simultaneous  congruences  has  been 
discussed  and  proved  in  reference  [6].  That  is; 

x  =  a,(mod(/;z,))  ;  /  =1 , 2,  3,  ...,  r.  (3-2-17) 

Then  to  have  a  solution  which  is  valid  for  any  pair,  one  has 

ai  =  aJ  mod  [g.c.dim^m,)) .  (3-2-18) 


(3-2-13) 

(3-2-14) 


This  results  in  a  single  solution  for  the  modulus 

Mr  =l.c.m  (ml ,  m2 ,  •  •  • ,  mr ) . 


(3-2-19) 
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For  example,  x  =  7mod(42)  and  x  =  15mod(5 1)  do  not  have  a  solution  since  g.c.d( 42,5 1)  = 
3  and  7  *  15mod(3). 

According  to  the  above  theorem,  there  is  a  unique  solution  to  these  congruences  for  a 
modulus  that  is  equal  to  the  product  of  all  the  given  ones.  The  first  known  source  of  such  a 
theorem  exists  in  the  Arithmetic  of  the  Chinese  writer  Sun-Tse  and  the  resulting  fonnula 
is  often  called  the  Chinese  remainder  theorem.  One  begins  by  forming  the  product 

M  =  mlm2  •  •  •  mr  (3-2-20) 

of  the  relative  prime  modulus  of  the  set  of  congruences.  When  M  is  divided  by  mt ,  the 
quotient 

M 

—  =  m,  •  •  •  m  (3-2-21) 

m] 

is  the  number  divisible  by  all  modulus  which  are  relatively  prime  to  m1 .  Similarly  the 
number 

M 

—  =  m]  •  •  •  mi  ,m+l  •  •  •  mr  (3-2-22) 

mi 

is  divisible  by  all  modulus  except  mi ,  to  which  it  is  relatively  prime.  For  each  i,  one  can 
solve  the  linear  congruence 


b,  —  =  1  mod(m  ) . 

v  ' 


(3-2-23) 


The  Chinese  remainder  theorem  can  be  stated  as:  Let  the  set  of  simultaneous  congruences 
given  for  the  modulus  mj  be  relative  primes.  Then 


x  =  at  mod(/u/)  ;  for  i  =1,  2,3,  ..., 


(3-2-24) 


For  each  i  one  determines  bi  through  the  linear  congruence 


bi  —  =  \  mod(/n( ) , 
m, 


Mr  =  l.c.rn  (ml ,  m2 ,  •  •  • ,  mr ) . 


(3-2-25) 

(3-2-26) 


The  solution  of  the  set  of  congruences  is  then 


x  = 


aA 


M 


m. 


■  +  a2b2 


M  M 

—  +---a  b  — 

nu  mr  j 


mod(  M ) . 


(3-2-27) 
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The  following  example  corresponds  to  the  three  congruences  with  x=2mod(3),  x  = 
3mod(5),  x  =  2mod(l).  If  M  =  105  and  then 


M_ 


=  35, 


M_ 

m3 


=  15. 


The  set  of  linear  congruences  will  be 

35bl  =lmod(3),  21  b2  =lmod(5),  25 b3  =lmod(7) 
and  it  has  the  following  solution 

bx  =  2 ,  b2  =  1 ,  b3  =  1 . 

Therefore 

jc  =  (3  -  2  -  35  +  3-1-21  +  2-1  •  15)mod(l05)  =  233mod(105)  =  23mod(105). 

The  Chinese  remainder  theorem  can  accurately  estimate  the  Doppler  frequency  and 
range  of  a  target  when  used  in  multiple  PRF  systems  as  compared  to  the  other  methods. 
In  addition  to  that,  other  methods  like  FFT  or  the  matrix  pencil  method  should  be  used 
before  applying  the  Chinese  remainder  theorem  to  estimate  the  Doppler  frequencies  for 
each  PRF.  The  problem  with  the  Chinese  remainder  theorem  approach  when  applied  to  a 
multiple  PRF  system  is  that  a  small  range  error  on  a  single  PRF  can  cause  a  large  error  in 
the  resolved  range  and  there  is  no  indication  that  this  has  happened.  Trunk  and  Brockett 
[5]  introduced  a  clustering  algorithm  to  resolve  the  range  and  the  Doppler  in  a  multiple 
PRF  systems. 
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Figure  3.2.3  How  Doppler  ambiguity  can  be  resolved  by  the  clustering  algorithm. 
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To  deal  with  this  problem,  one  should  be  given  a  measure  of  the  error  at  the  estimated 
values  of  the  Doppler.  As  shown  in  Figure  3.2.3,  the  Clustering  algorithm  calculates  the 
error,  denoted  by  C(  j) ,  between  the  estimates  and  the  average  of  the  estimates  of  Doppler 
that  is  the  same  as  the  calculation  of  variance.  At  the  actual  Doppler  frequency,  the 
variance  by  the  PRFs  should  be  the  minimum.  Without  any  noise,  the  C(j)  should  be  zero 
at  the  actual  Doppler  frequency.  The  value  of  C(j)  at  the  original  Doppler  will  also  be  the 
same  as  the  variance  for  the  noise  in  the  measurement  and  the  maximum  may  be  obtained 
by  the  Cramer-Rao  bound.  The  average  squared  error  C(  j)  for  m  consecutive  Doppler  is 

1  J+m ,  _  ,  9 

c(,/)=- 

m  i=j+l  1 

(3-2-28) 

where  /  =  the  average  value  of  the  m  ordered  Doppler, 
m  =  number  of  the  PRF, 

foi  =  tested  frequency  for  z-th  PRF  and  o  times  the  folded  one. 

Consider  a  Doppler  frequency  of  25.5kHz  to  illustrate  how  the  algorithm  works  for  the 
multiple  PRF  systems.  If  the  three  PRFs  are  PRFl  =  5kHz,  PRF ,  =  6kHz,  PRFi  =  7kHz, 

then  C(j)  can  be  calculated  and  is  shown  in  Figure  3.2.4.  C(j)  is  minimized  at  the  target 
frequency  of  25.5kHz  which  is  equal  to  the  argument  j  =  63. 


Figure  3.2.4.  The  result  of  applying  the  clustering  algorithm. 
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3.3  LEAST  SQUARES  METHOD 


The  concept  of  spectral  analysis  to  nonunifonnly  sampled  data  using  Least  squares  was 
first  proposed  by  Vanicek  in  1970  [10].  Lomb  (1975)  developed  this  method  and  showed 
that  a  correlation  exists  between  the  height  of  the  spectrum  at  any  two  frequencies  which 
is  equal  to  the  mean  height  of  the  spectrum  due  to  a  sinusoidal  signal  of  frequencies  f 
and  f2  [11].  These  correlations  reduce  the  distortion  in  the  spectrum  of  a  signal  affected 
by  noise  which  is  an  additional  benefit  to  using  unevenly  spaced  data  [11].  Further  studies 
have  been  done  by  Scargle  (1982)  in  which  he  provided  a  simple  estimate  of  the 
significance  of  the  height  of  a  peak  in  the  power  spectrum  through  the  false  alarm 
probability  [19].  Feraz-Mello  used  non-orthogonality  of  the  basis  functions  when  the 
sampling  is  uneven  and  then  applied  the  Gram-Schmidt  orthogonalization  procedure 
which  is  basically  equivalent  to  a  periodogram  based  method  [21]. 

The  periodogram  approach  to  the  evaluation  of  the  spectrum  fonn  a  set  of 
nonunifonnly  sampled  data  provides  a  scan  of  a  given  frequency  range.  This  is  obtained 
by  fitting  sines  and  cosines  functions  in  a  Least  squares  fashion  to  the  data  and  plotting  the 
correlation  of  the  data  for  each  frequency.  The  Least  square  spectrum  provides  the  best 
measure  of  the  power  contributed  by  the  different  frequencies  to  the  overall  variance  of 
the  data.  Therefore  this  can  be  regarded  as  the  natural  extension  of  the  Fourier  methods  to 
nonunifonnly  spaced  data.  The  frequency  increment  can  be  determined  with  any  precision 
and  that  is  an  additional  benefit  of  using  this  method.  Additional  advantages  can  be 
derived  from  the  uneven  or  random  sampling  which  is  absence  of  aliasing  if  the  sampling 
were  to  be  completely  random  [9,  13,  15].  It  is  known  that  in  such  situations  the  spectrum 
would  be  completely  alias  free  [12]. 

Even  though  the  periodogram  analysis  has  many  benefits,  it  also  has  some  drawbacks. 
For  example  it  cannot  evaluate  the  spectrum  for  negative  frequencies  since  it  is  a  power 
spectrum  of  a  real  sequence.  The  estimated  peak  does  not  precisely  correspond  to  the  true 
magnitude  of  the  signal.  The  error  in  the  peak  is  mainly  from  the  nonuniform  spacing, 
which  does  not  have  the  same  source  of  error  as  in  the  FFT  in  which  the  error  is  primarily 
due  to  the  finiteness  of  the  sequence.  A  fonnulation  that  can  resolve  both  positive  and 
negative  frequencies  without  loosing  any  of  the  benefits  of  the  periodogram  approach  has 
been  studied  in  this  section.  Additional  properties  of  the  Least  squares  approach  have  been 
investigated  in  section  3.3.2  to  reduce  the  number  of  computations  in  a  real  time 
operation. 

3.3.1  Formulation  of  the  Least  Squares  Method 

Let  a  continuous  complex  signal  x{t)  be  sampled  at  time  instants,  t  =  tk ,  for  k  =  0,  1,2, 

. . .,  TV— 1 .  We  are  interested  in  looking  for  a  harmonic  component  of  frequency  co ,  so  that 

h(t)  =  {a  +  ib)cos[co(t  -  r)]  +  (c  +  id)sin[a)(t  -  r)]  (3-3-1) 

where  a,  b,  c,  d  and  t  are  real  constants.  The  delay  parameter  t  enables  one  to  select  any 
arbitrary  location  of  the  origin  in  the  time  axis.  To  estimate  a,  b,  c  and  d  the  following 
mean  square  difference  is  minimized, 
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F  =  ^  |  x(tk )  -  {a  +  ib)  cos[®(ft  -  r)]  -  (c  +  id)  sin[<z>(tA,  -  r)]  |  (3-3-2) 

k=0 

with  respect  to  the  unknowns.  Taking  derivative  of  F  with  respect  to  the  unknowns  a,  b,  c 
and  d  will  produce  the  normal  equations  which  are  of  the  form 


^  =  0 
da 

or  equivalently 

7V-1  N- 1 

EJ(0  cosMa  -r)]+Xx(A)  cosK4-^)] 

k=0 

N-l  N-\ 

=2aT,  cos2  [a>(tk -r)]+2c^cos[&>(4 -r)]  sin [o)(tk-r)]9  (3-3-3) 


k=0 


k=0 


k=0 


where  x(tk)  is  the  complex  conjugate  of  x(tk).  Since  x  is  a  free  parameter,  it  is  selected 
so  as  to  simplify  the  normal  equations;  that  is, 


N-l 


^cos[&>(4  -  r)]sin[^(^  -  r)]  =  0 . 


k=0 


Solving  for  x  will  yield 


tan(2cor)  = 


N-l  /  \ 

X  sm(2a)  fJ 


X  cos(2o%) 
k  =  0 


Similarly  for  the  parameters  b,  c  and  d ,  we  enforce 

dF  dF  dF 


db  dc  dd 


=  0. 


N-l 


k=0 


N-l 


k= 0 


N-l 


k= 0 


N-l 


k=0 


k= 0 


k= 0 


(3-3-4) 


(3-3-5) 


Use  of  (3-3-4)  will  yield  the  following  equations 

N-l  N-l 

^/x^Jcos^.  -rjj-^zx^Jcos^^  -r)]  =  2h^cos2[ru(tA.  -r)],  (3-3-6) 


Z*(OsinW*  -r)]  +  Zx(^)sinW^  -r)l  =  2c^sin2 [co(tk -t)],  (3-3-7) 


26 


(3-3-8) 


N- 1 


X4a  )sinh(^  -  r)] -  X'4  )sinh(^  -  T)]  =  2^Xsin2  [®fat  -  r)]  • 


k=0 


k= o 


A:=0 


The  resulting  values  are: 


N- 1 


W-1 


(2  = 


X  Re[-4 )]  cos[4a  -  *■)]  X  ImIx(^ )]  cosb(^  -  t)\ 

- ,  b= k=0 


k=0 


N-l 


X  cos 2  [4a  -r)] 


N- 1 


X  cos2  [4a  -r)] 


&=0 


£=0 


JV— 1 


7V-1 


C  = 


k= 0 


AM 


and 


X  Re[4 )]  ~t)]  X  Im[4 )] sui[4a  -  t )] 

^  =  ^4=1 - > 

X  sin 2  [4a  ~T)\ 

k= 0 

AM 

X^Jsinkf,  ~t)] 


Xsin2[^  -r)] 


k= 0 


AM 


Xx(Oc°s[4t-T)] 


a  +  ib  =  *4 


N- 1 


c  +  id  = 


k= o 


Xc°s2[4-  -r)] 


AM 


X  sin 2  [4a-  -0] 


*=0 


&=0 


Substituting  (3-3-10)  into  (3-3-1)  will  yield 


AM 


N- 1 


(3-3-9) 


(3-3-10) 


X4a  )cos[4*  -  r)]  X4  )sin[®h  -  r)] 

4)  =  - cos[4  -  r)]  +  2=4 - sin[4  -  r)]  (3-3-1 1) 


X  cos2  [4a  -01 


X  sin 2  [4a  -0] 


k=0 


k= 0 


The  power  in  the  harmonic  component  at  frequency  co  is  given  by 


44X1  4 

k= 0 


.  N- 1 


.  N- 1 


+ 


/b|  X c°s2  [4a  ~  r)]  + 1  c  +  id\  ^ sin2 [&>(^  - t)\ 


k=0 


k= 0 


AM  |  |  iv-i 

E  4a ) cos[4*  -  r)j  \  \  X  x(tk ) sin[4*  -  r 

A=0 


Zcos2[4a  -0] 

A=0 


-  + 


JV— 1 

E; 

k-0 


Esin2[4/f  -r)] 

A=0 


(3-3-12) 


Equation  (3-3-12)  is  a  complex  version  of  the  Lomb  periodogram  [11].  Observe  that  (3-3- 
12)  yields  the  same  value  for  co  and  -co  since  it  is  squared.  Obviously,  this  expression  is 
not  suitable  for  negative  frequencies.  To  obtain  the  phase  component  of  the  spectrum  from 
the  power  representation  (3-3-12),  a  frequency  response  E(co  )  of  the  following  form  is 
assumed: 
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(3-3-13) 


N-l  N-l  - 

£W  =  Z  a  cos[&>  (tk  -  t)  +  / 3  sin[<y  (tk  -  r)  , 


k= 0 


1=0 


where  r  is  a  free  parameter  as  defined  in  (3-3-1).  The  corresponding  power  spectrum  can 
be  written  as 


i  ,  „  ,  9 

r  , 

2 

v^1  r  -I 

= 

^]acos|ft>  (t/c  -  r)J 

+ 

Z^sinh  (fc  -  *■)] 

k=0 

k= 0 

(3-3-14) 


Matching  (3-3-12)  and  (3-3-14)  for  all  x(tk)  will  give  the  unknown  coefficients  a  and  [j . 


i.e.. 


Thus, 


N-l  , 

^  a  cos  co  (tk  -  r)J  =  ± 


N-l  , 

^x(4)cos[m  (tk-T) \ 


k= 0 


k= 0 


N-l 


N-l 


Z^sin[®  (t,-r)]=  + 


Zcos2[m  (tk  -  r)] 

V  k= 0 

N-l  r  , 

Ex(4)sin[®  (tk-r) \ 


k=0 


k= 0 


Zsin 2[v  (f*-r)] 


k=0 


N-l 


N-l 


N-l  N-l 

Zc°s2  i®  (h  -  *■)]  i  Zsin>  (**  ~  r)J 


k= 0 


i=0 


(3-3-15) 


2>(0  cos  [o){tk-T)\  sin[ru  (tk  —  r )] 

E(cq)  =  ±^='  - ±i-*%=  =-  ■  (3-3-16) 


There  can  be  four  possibilities  of  choosing  a  sign  and  the  problem  is  how  to  choose  the 
correct  one.  The  expression  for  the  correct  frequency  can  easily  be  obtained  from  the 
analogy  of  the  conventional  Fourier  transfonnation  or  by  inserting  the  following  test 
signal  x(tk)~  ei0,,h  and  observing  (3-3-16)  at  the  frequency  cox  and  ~(ol.  The  resulting 
expression  takes  positive  signs  for  both  the  terms,  i.e., 

N-l  j.  1  :Y-|  p  , 

Z )cosl.^  -  r)J  Z  )sink  (**  -  T)\ 

e{co)  =  ^=  —  +  j*=(j— _  (3-3-17) 

J  Z cos2  [®  ({k  - *■)]  J  Z sin2  [®  (lk  -  t)] 

V k= 0  V  k= 0 

The  phase  information  can  be  obtained  using  this  expression. 
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To  illustrate  the  significance  of  (3-3-12)  and  (3-3-17),  consider  a  signal  of  the  form 

/( xk )  =  e2'27ttk'  +  2e~2SlMti  -  2.5e*2Mki , 

where  tk is  a  random  number  uniformly  distributed  between  [0,  100]  for  k  =  1,2,  ...  ,  100. 
Observe  that  this  is  the  same  signal  as  in  the  section  3.1  except  for  the  average  sampling 
frequency  is  much  higher  than  the  previous  one.  Now  (3-3-12)  and  (3-3-17)  are  used  to 
estimate  the  spectrum  of  the  nonuniformly  sampled  data  and  the  result  is  shown  in  Figure 
3.3.1.  Average  sampling  frequency  is  1  FIz  which  is  much  lower  than  the  Nyquist 
sampling  rate  (4  Hz  since  this  is  a  complex  signal).  Here,  the  sampling  frequency 
corresponds  to  the  highest  vale  of  the  signal  and  not  twice  the  highest  frequency  as 
required  for  the  conventional  Nyquist  sampling  theorem  to  hold.  As  seen  in  Figure 
3.3.1(b)  the  negative  frequency  component  located  at  -3.5  Hz  is  now  distinguishable 
through  the  use  of  Equation  (3-3-17).  This  is  not  possible  in  Figure  3.3.1(a)  which  is  the 
original  formulation  for  the  Lomb  periodogram.  Note  that  the  absolute  value  has  been 
taken  and  squared  in  Figure  3.3.1(b). 
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Figure  3.3.1  Comparison  between  the  Lomb  periodogram  and  the  new  modified  scheme. 

3.3.2  Hilbert  Transform  Relationship 

We  can  reduce  the  computation  time  in  the  evaluation  of  the  spectrum  time  by  using  a 
Hilbert  transfonn  relationship  between  the  real  and  the  imaginary  parts  of  (3-3-17),  E{co)- 
This  means  that  the  time  domain  response  of  E(co)  is  causal  and  this  is  shown  next.  The 
Hilbert  transformation  can  be  obtained  by  performing  two  Fast  Fourier  transforms  and  one 
multiplication.  Assume  that  the  number  of  frequency  steps  is  M,  and  then  the  operation 
count  for  the  Hilbert  transform  will  be  2M\ogM + M .  The  operation  count  for  the 


(b)  Modified  one  using  (3-3-17) 


~i - 1 - 1 - r~ 


~i - 1 - 1 - r 


(a)  Lomb  Periodogram  using  (3-3-12) 
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evaluation  of  Equation  (3-3-17)  is  4{N+ 1  )M  where  N  is  the  number  of  time  domain  data 
samples.  By  utilizing  the  Hilbert  transform  relationship  the  real  part  of  (3-3-17)  can  be 
obtained  from  the  imaginary  part,  and  vice  versa  with  reasonable  accuracy  and  the 


processing  time  will  be  reduced  by  a  factor  of 


2  log(M)+ 1  +  (2N  +  2) 
4N  +  4 


approximately. 


If 


M  and  N  are  large  numbers,  a  maximum  of  50%  of  reduction  in  the  computation  time  by 
utilizing  the  Hilbert  transformation  is  obtained. 

Assume  a  causal  time  domain  signal  x{t)  which  exists  on  [0,«]  where  a  is  a  finite 
number,  then 


EW  = 


(3-3-18) 


Here,  the  tenn  v  is  ignored  since  a  is  assumed  to  be  a  large  number  compared  to  the 
period  of  the  signal  and  it  is  assumed  that  x{t)  is  a  time  invariant  sequence.  When  the  time 
step  is  small,  we  can  replace  the  summation  in  (3-3-18)  by  an  integration,  i.e., 


j"x(t)cos(<u  t)  dt  Jx(t)sin(<u  t)  dt 


(3-3-19) 


Application  of 


Jeos2 

0 


( 'co 


sin  (2  coa) 
4co 


,  J sin2(<w  f)  dt 

o 


a  sin(2  coa) 
4 a 


will  transfonn  equation  (3-3-19)  to 


E(co)  = 


a  a 

Jx(t)cos(<z>  t)  dt  |x(t)sin(<z)  t)  dt 


ja  ~  sin(2 coa) 
2  4  co 


+  l 


;  0 


fa  sin(2  coa) 
2  4  co 


(3-3-20) 


If  we  assume  \a>a\  » 1 ,  then  sin(2  coa)  «  j2coaj  and 


/2,  a  2  a 

— |x(t)cos(ry  t)  dt  +  ij—j x(t)sin(co  t)  dt .  (3-3-21) 
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Since  x(t)  is  causal,  the  real  and  imaginary  part  of  E(co)  will  be  related  by  the  Hilbert 
transform  relationship  when  |<wa|»l.  Figure  3.3.2  compares  Ex  {co)  with  the  Hilbert 
transform  of  E2{co)  where  E{co)  =  E]  {co)  +  iE2 {co)  and  E{{co),  E2{co)  are  real.  Figure 
3.3.2(a)  corresponds  to  the  case,  when  \a>a\  is  relatively  a  large  number,  and  then  the  two 
results  will  coincide  with  each  other.  When  \coa\  is  small,  Figure  3.3.2(b)  shows  that 

there  are  some  differences  between  the  two  curves  where  co  is  of  a  small  value.  To 
illustrate  the  applicability  of  this  method  we  consider  the  same  signal  as  in  the  previous 
example: 


f(xk )  =  e2'2Mti  +  2e~3S2Mki  -  2.5e^ki  ;  0  <tk<a 

and  the  Matlab  function  HILBERT  is  used  to  compute  the  Hilbert  transform  of  E2  {co) . 
Processing  time  to  obtain  the  estimator  (3-3-17)  can  be  measured  by  changing  the  number 
of  frequency  steps  and  the  number  of  time  domain  data  samples.  The  result  is  shown  in 
Figure  3.3.3.  By  utilizing  the  Hilbert  transformation,  the  processing  time  has  been 
reduced  by  46%. 


(a)  When  a  is  large  (a  =  2) 


(a)  When  a  is  small  (a  =  0.5) 


Frequency  [Hz] 


Figure  3.3.2  El{oo)  and  Hilbert  transform  of  E2{co)  for  different  value  of  (OCX  .  If  (OCX  »1  results 

coincide  with  each  other  as  shown  in  (a). 

3.3.3  Estimation  of  the  amplitude 

As  seen  in  Figure  3.3.1,  none  of  the  Lomb  periodogram  methods  or  the  modified  ones 
provide  the  exact  amplitude  of  the  signal.  The  error  obviously  comes  from  the  uneven 
spacing  and  the  aliasing  between  the  different  frequencies. 
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Since  the  estimates  of  the  frequencies  giving  rise  to  the  peaks  are  not  much  different 
than  the  actual  ones,  we  can  estimate  their  amplitudes  from  these  frequencies  by  using  a 
Least  square  method.  If  the  signal  is  a  sum  of  exponentials,  then 

4)=iv“"^=  1,2,  ..., N,  (3-3-22) 

/= 1 


where  col 

Ah) 

L 

and  Al 


frequencies  which  give  highest  peaks, 

given  data  with  respect  to  unevenly  spaced  point  tk  , 

number  of  frequency  components, 

unknown  magnitudes. 


CPU  time  (a)  CPU  time  reduces  by  46%  when  using  the  Hilbert  transform 


CPU  time  (b)  CPU  time  reduces  by  47%  when  using  the  Hilbert  transform 


Figure  3.3.3  Processing  time  is  reduced  by  using  the  Hilbert  transform  relationship. 


By  rewriting  (3-3-22)  as 
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(3-3-23) 


A]) 

eic°\h 

ei0}2h  . 

1 

•§ 

Ai) 

ei»\h 

eic°ih  . 

.  eimLh 

An). 

e'®ilv 

eia>2tN  _ 

.  g  iO>LfN 

A 

A 

A 


and  using  the  pseudo  inverse,  a  vector  A  containing  all  the  amplitudes  can  be  obtained 
from  the  unevenly  sampled  points  of  the  signal  x(tk)  corresponding  to  the  estimated 


frequency  cok  as, 


where  B  = 


ico2t\ 


io)xt  2 


io)2t2 


A  =  (b‘b)'  B*  f 


(3-3-24) 


eir»L‘\ 

2 


and  B  is  a  conjugate  transpose  of  B. 


gic0l  *N  giC°2  fN 


NxL 


The  same  signal  as  described  by  the  first  example  has  been  used  to  verify  (3-3-24)  and  the 
result  is  shown  in  Figure  3.3.4.  The  three  signal  components  with  unit  amplitudes  can  be 
obtained  precisely  by  utilizing  (3-3-24)  while  the  amplitudes  obtained  from  using  (3-2-17) 
have  some  differences  when  compared  with  the  true  value. 
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Figure  3.3.4  Result  of  Equation  (3-3-19).  Since  the  periodogram  does  not  give  accurate  values  of  the 

amplitudes,  (3-3-19)  can  be  used. 
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3.3.4  Summary 

Unevenly  spaced  spectrum  using  the  Least  squares  method  has  been  studied  in  this 
chapter.  The  well-known  Lomb  periodogram  approach  has  many  benefits  but  it  cannot 
discriminate  between  positive  and  negative  frequencies.  Using  a  modified  scheme, 
positive  and  negative  frequencies  can  be  discerned  without  losing  any  of  the  benefits  of 
the  Lomb  periodogram.  The  method  to  estimate  the  magnitude  of  the  signal  using  the 
periodogram  approach  has  also  been  described. 

One  of  the  properties  of  the  periodogram  approach  is  that  a  Hilbert  transform  pair 
relates  the  real  and  imaginary  parts  of  the  coefficients.  By  utilizing  this  property,  the 
computation  time  for  the  spectrum  can  be  reduced  by  half. 

3.4  MULTI-RESOLUTION  ANALYSIS 

In  an  effort  to  reconstruct  a  band-limited  signal  from  an  unevenly  sampled  signal  using 
multirate  filter  banks,  a  multi-resolution  approach  has  been  taken  by  Vaidyanathan  [25,  26 
and  27].  ft  is  possible  to  reconstruct  a  band-limited  sequence  from  a  data  sequence 
followed  by  decimation. 

Since  the  spacing  of  the  multiple  PRF  is  not  uniform  and  is  not  totally  random,  it  can 
be  characterized  by  an  unevenly  decimated  version  of  the  signal.  In  that  case,  one  can 
utilize  quadrature  mirror  filters  (QMF)  to  get  the  original  signal.  Two  preliminary 
conditions  should  be  satisfied  when  using  this  approach  to  reconstruct  a  signal  from  an 
unevenly  sampled  data.  First,  the  signal  should  be  band-limited.  In  real  life,  there  cannot 
be  a  band-limited  sequence  since  our  sample  length  in  time  is  finite.  If  numerous  data 
points  are  obtained  or  the  signal  spectrum  is  concentrated  in  a  certain  region,  then  one  can 
consider  the  signal  to  be  approximately  band-limited  for  all  practical  purposes.  Second, 
the  sampling  frequency  cannot  be  totally  random  but  it  should  have  a  common  interval  so 
that  each  sample  in  the  multiple  PRF  signal  is  an  integer  multiple  of  the  common  interval. 
Since  multiple  PRF  system  satisfies  the  second  condition,  the  main  contribution  to  the 
error  would  be  the  finiteness  of  the  data  samples  and  the  accuracy  in  the  design  of  the 
filters. 

This  method  illustrates  the  correlation  between  sampling  rate  and  the  bandwidth.  If  the 
bandwidth  becomes  large,  the  number  of  samples  should  increase  within  a  given  period  of 
the  signal  so  as  not  to  loose  any  information.  Conversely,  if  the  number  of  data  samples 
becomes  small,  like  missing  data  or  due  to  decimation  of  the  data,  the  bandwidth  should 
decrease  proportionally.  In  an  ideal  case,  the  QMF  analysis  can  recover  the  unknown 
samples  in  between  the  missing  data  exactly.  A  two  PRF  case  has  been  presented  in  this 
section  and  this  can  be  extended  to  the  multiple  PRF  case. 

3.4.1  Two  PRFs  Case  (20kHz  and  30kHz) 

It  is  known  that  if  a  sequence  is  band-limited  and  the  average  sampling  rate  satisfies 
the  Nyquist  sampling  criterion,  then  the  sequence  can  fully  be  recovered  regardless  of  the 
sampling  rate  whether  it  be  unifonn  or  not.  Assume  that  we  have  M  uniform  samples  in  a 
period,  then  the  maximum  angular  frequency  content  of  the  signal  would  be  \co\  <  n  .  If 
we  have  a  decimated  sequence  which  has  N  randomly  chosen  samples  out  of  M  uniform 
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I  I  N 

samples,  then  the  bandwidth  of  the  sequence  should  be  reduced  to  \co\<n —  for  a 

M 


complete  recovery  of  the  signal  without  aliasing. 

Figure  3.4.1  provides  an  example  of  the  data  sampled  by  two  PRFs,  20kHz  and  30kHz. 
Since  the  greatest  common  divisor  of  those  two  frequencies  are  60kHz,  the  sequence  is 
repeated  every  6  intervals  (every  4th  sample).  We  want  to  recover  the  missing  2  samples 

i  i  4 

occurring  at  time  step  1  and  5  out  of  6  samples  from  a  band  limited  sequence  of  \co  <  —  n  . 

6 

Perfect  reconstruction  can  be  achieved  by  a  QMF  system  shown  in  Figure  3.4.2  if  we  can 
derive  a  set  of  filter  banks  which  can  carry  out  a  perfect  reconstruction  of  X(z) .  Hence, 
the  missing  samples  can  be  recovered.  The  intermediate  step  labeled  U t  s  in  Figure  3.4.2  is 


the  unevenly  sampled  sequence.  One  can  estimate  the  output  x(z)  using  Equations  (3-4- 
1)  and  (3-4-2).  The  output  z  transformation  of  the  M-fold  decimator  (up  sampler)  Y(z)  is 
[25]: 


Y(z) 


M 


(3-4-1) 


Figure  3.4.1  Data  generated  by  sampling  a  signal  using  2  PRFs  (20kHz  and  30kHz). 


Figure  3.4.2.  Block  diagram  for  the  reconstruction  of  the  band-limited  signal  when  sampled  by  the  2  PRFs 

(20kHz  and  30kHz). 
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The  output  Z  transform  of  the  Z-fold  expander  (down  sampler)  Y(z)  is  given  by  [25]: 


_jln /  _jn/ 

We  define  W  =  e  'M=e  73 
equation  (3-4-1),  as 


Y{z)  =  x(zl). 


(3-4-2) 


and  the  Ut  in  the  Figure  3.4.2  would  be  obtained  using 


6  k=0  V  ’ 

u&)  =  \p*y*wk)  2*(z%Wkj,  (3-4-3) 

U3(z)  =  l'E(z^wk)  3^z^Wkj, 

i/4(z)  =  ^(ZV) 

By  using  (3-4-2),  the  Vs  are  given  by 

*i(4  =  7 

6£=0 

F2(z)  =  i  UzWkY x{zWk)H2,  (3-4-4) 

o  ^=0 

r3W=7  l:(z»ft)r3x(zrt)//3, 

6/t=0 

F4(z)  =  ^  l(z^)'V(z^)//4. 

6/t= 0 

The  output  X(z)  will  be  the  sum  of  all  the  Vi  ’s  and  this  should  be  a  purely  delayed 
version  of  the  original  sequence  if  it  were  to  be  a  perfect  reconstruction  of  our  signal. 
Thus, 

X(z)=VI(z)+V2(z)+V,(z)+Vt(z) 
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1 

—  < 

6 


\[x(z)  +  X(zW)  +  x(zW2)+  x(zW3 )+  x(zW4 )+  x(zW 5 )]//, 

X(z)  +  W-2X(zW)  +  W-4x(zW2)+W-6x(zWi)+W-*x(zW4)+W-l0x(zW5 )]z  2H2 
X(z)  +  W-3X(zW)+W-6x(zW2)+W-9x(zW3)+W->2x(zW4)+W-15x(zW5)]z-3H3 
.X{z)  +  W-4X{zW)+W-8x(zW2)+W-nx(zW3)+W-l6x(zW4)+W-20x(zJV5  )]z~4/74 


=  z-kX(z),  (3-4-5) 


where  k  is  any  constant.  The  aliased  terms,  X(zW) ,  x(zW2) ,  x{zW3) ,  x{zW4)  and 
x(zW') ,  should  be  zero  for  perfect  reconstruction  and  the  output  signal  term  should  be  a 
delayed  version  of  the  original  signal.  Matching  term  by  term  in  each  of  the  two 
expressions  will  yield 

—X(z)  [Hl  +z  2H,+z  3H3  +  z~4H4\  =  z~k X(z), 

6 

(3-4-6-a) 

-X(zW)  [hi+W  2z  2H2  +  W  3z  3H3  +W~4z~4H4]=  0, 

6 

(3-4-6-b) 

-x(zW2)  [hx+W  4 z~2H2  +  W  6z-3H3+W  8z-4H4}=0, 

6 

(3-4-6-c) 

^ x(zW 3)  [h,  +W~6z-2H2  +W~9z~3H3+W-nz-4H4\=0, 

(3-4-6-d) 

-x(zW4)  [hx  +  W~*z-2H2+W-12z-3H3  +W~16z-4H4\=  0, 

6 

(3-4-6-e) 

-x(zw5)  [hx  +w-10z-2h2  +w-'5z-3h3  +w-20z-4h4]=o  . 

(3-4-6-f) 

We  have  4  unknowns,  Hl,H2,H3  and  H4,  and  6  equations  which  cannot  be  solved 
simultaneously  in  general.  In  our  case,  the  shifted  terms  X{zW) ,  x{zW2) ,  x{zW3) , 
x[zW 4)  and  x(zW5)  exist  only  in  limited  regions  in  the  frequency  domain  and  we  can 
make  those  equations  satisfy  the  constraints  in  each  region.  Figure  3.4.3  shows  the  shifted 
versions  of  the  sequence  X(z) .  In  region  1  of  Figure  3.4.3,  Equation  (3-4-6-a),  (3-4-6-b), 

(3-4-6-c)  and  (3-4-6-d)  need  to  be  satisfied  since  x(zW4)  and  x(zW5)  are  zeros  and  it 
satisfies  all  6  equations  in  region  1 .  Solving  the  equations  will  give  the  filter 
characteristics  in  region  1.  The  filters,  Hx ,  z  2 H2 ,  z  3 II 3  and  z  4/74 ,  should  be  constant 
numbers  for  each  frequency  region  since  W  is  a  constant.  In  region  2,  Equation  (3-4-6-a), 
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(3-4-6-b),  (3-4-6-c)  and  (3-4-6-f)  need  to  be  satisfied  since  x{zW 3)  and  x{zW4)  are 

zero  and  it  satisfies  all  the  six  equations.  Solving  the  equations  will  give  the  filters 
characteristics  in  the  region  2.  In  region  3,  Equation  (3-4-6-a),  (3-4-6-b),  (3-4-6-e)  and  (3- 

4-6-f)  need  to  be  satisfied.  x(zW2)  and  x(zW3)  is  zero  which  makes  it  possible  to 

satisfy  all  the  six  equations,  simultaneously.  Solving  the  equations  will  give  the  filter 
characteristics  in  the  region  3.  In  region  4,  Equation  (3-4-6-a),  (3-4-6-d),  (3-4-6-e)  and  (3- 
4-6-f)  need  to  be  satisfied  since  X(zW)  and  x{zW2)  are  zero  and  this  satisfies  all  the  6 

equations.  Solving  all  these  equations  will  give  the  filter  characteristics  in  the  region  4. 

Once  the  filter  coefficients  for  each  region  are  calculated,  they  should  be  combined  to 
make  a  set  of  filters.  Note  that  each  filter  from  the  equation  yields  Hx,  z  2H2,  z  77, 

and  z  4H4  and  the  delay  should  be  taken  out  from  them  when  realizing  the  filters.  The 
results  of  those  filters  are  shown  in  Figure  3.4.4. 


(a)  The  spectrum  of  a  signal 


(b)  Shifted  versions  of  the  spectrum 


Frequency  [rad] 


Figure  3.4.3  Spectrum  of  a  sampled  signal  and  its  shifted  versions. 
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Real  part  of  Hi 


Real  part  of  z"2H2 


Imaginary  part  of  Hi 


Real  part  of  z'3H3 


Real  part  of  z"4H4 


Imaginary  part  of  z"3H3 


Figure  3.4.4  Filters  to  generate  a  perfect  reconstruction  of  the  signal  for  the  system  of  Figure  3.4.2. 


A  numerical  example  outlined  in  Figure  3.4.2  utilizing  the  system  of  filters  as 
described  in  Figure  3.4.4  is  used  to  illustrate  how  this  method  performs.  The  original 
signal  is  described  by 

f{xk)  =  e2'27a>l  +2e^5'lnXkl  -2.5e4'2nXki ;  for  k  =  1,  2,  ...,  168.  (3-4-7) 

The  signal  is  sampled  at  5.5556FIz  and  8.3333FIz  (equivalent  sampling  at  20kFIz  and 
30kHz  with  time  scaling)  simultaneously  as  shown  in  Figure  3.4.1  and  it  is  equivalent  to 
sampling  the  signal  at  16.6667  (60kHz)  and  then  unevenly  decimating  it.  Note  that  the 
maximum  frequency  in  the  signal  is  4Hz  which  is  less  than  half  the  equivalent  sampling 
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frequency  5. 5556Hz  ( =N/M  •  16.6667/2).  The  frequency  domain  result  is  given  in  Figure 
3.4.5.  Observe  that  the  original  signal  is  not  strictly  band-limited  due  to  the  finite  number 
of  data  samples  but  the  transformed  one  is  band-limited.  Figure  3.4.6  produces  the 
resulting  time  domain  sequence. 

3.4.2  Multiple  PRF  Case 

The  same  procedure  as  described  for  the  2  PRF  case  can  be  applied  to  the  multiple  PRF 
case.  Define  M  as  the  total  number  of  time  steps  in  one  period  where  decimation 
repeatedly  occurs.  If  the  PRFs  are  relatively  prime  numbers  of  each  other,  the  time  step  M 
will  be  given  by  the  least  common  multiplier  ( l.c.m )  of  the  PRFs.  For  example,  the 
sampled  data  for  the  3  PRFs  3Hz,  4Hz  and  5Hz  which  are  relatively  prime  numbers  is 
repeated  every  60  time  steps  as  illustrated  in  Figure  3.4.7.  If  the  PRFs  are  not  relative 
prime  numbers,  then  M  should  be  divided  by  the  greatest  common  divisor  ( g.c.d)  of  the 
PRFs.  Sampling  the  data  by  PRFs  of  3Hz,  4Hz  and  5Hz  is  same  as  sampling  the  data  at 
30kHz,  40kHz  and  50kHz  except  that  the  scale  factor  and  M  can  be  recalculated  by 
dividing  the  PRFs  with  their  g.c.d,  i.e., 

M  =  l.c.m(PRFl ,  PRF2  •  •  •)  (3-4-7) 

g.c.d {PRFl ,  PRF2  •  •  •) 


400 
300 
200 
100 
0 

-3-2-1  0  1  2  3 

(b)  Spectrum  of  the  estimated  sequence 
400 

300 

200 

100 

0 

-3-2-1  0  1  2  3 

Frequency  [rad] 

Figure  3.4.5  Spectrum  (magnitude)  for  the  data  of  Figure  3.4.2  when  using  the  QMF  filters  as  shown  in 

Figure  3.4.4. 


(a)  Spectrum  of  the  original  sequence 
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(a)  Original  sequence 


0.4  0.6  0.8  1  1.2  1.4  1.6  1.8  2 

Time 

Figure  3.4.6  Time  domain  data  for  the  example  of  Figure  3.4.2  when  using  the  QMF  filters  as  shown  in 
Figure  3.4.4. 


The  number  of  known  samples  in  one  period  adds  up  to  the  PRFs  reduced  by  1  if  they 
are  relative  prime  numbers.  If  they  are  not  relative  prime  numbers,  it  can  simply  be 
divided  by  their  g.c.d.  as  before.  Define  N  as  the  number  of  known  samples  in  the  M  time 
intervals,  where 


'EPRFi 

N  = - - - 1 

g.c.d(PRFl ,  PRF2  ■  ■  •) 


(3-4-8) 


For  the  30kHz,  40kHz  and  50kHz  case  we  will  have  1 1  known  samples  in  one  period.  We 
have  49  (=  60  -  11)  unknown  samples  in  one  period  in  this  case  (represented  by  the  dots 
in  Figure  3.4.7).  The  ratio  of  the  number  of  unknown  samples  to  the  known  samples 
determines  the  maximum  bandwidth  of  the  signal  with  respect  to  the  Nyquist  rate.  The 
real  data  sampled  by  the  exact  Nyquist  frequency  or  at  a  higher  frequency  can  have  a 
maximum  bandwidth  of  2  times  the  Nyquist  frequency  for  an  exact  reconstruction.  A 
more  sparsely  sampled  data  than  the  Nyquist  rate  should  have  less  bandwidth  for  an  exact 
reconstruction.  The  maximum  bandwidth  allowed  for  the  signal  would  be 

i  i  N 

\eo\<x—.  (3-4-9) 

M 

Obviously,  the  bandwidth  if  sampled  by  the  Nyquist  rate  is  \co\  <  n . 
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o  =  30kHz 
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Figure  3.4.7  Sampled  data  generated  by  the  three  PRFs. 


Even  though  the  signal  is  band-limited  by  Equation  (3-4-9),  the  effective  sampling 
frequency  increases  to  the  l.c.m  of  those  three  PRFs  of  3,  4  and  5Hz  which  is  60Hz.  The 
actual  maximum  resolvable  Doppler  frequency  also  increases. 


The  number  of  filters  that  should  be  designed  is  the  same  as  the  number  of  known 
samples  (N)  as  in  the  Figure  3.4.2.  The  spectra  should  be  divided  into  M  regions  to 

2  k 

calculate  each  of  the  filter  coefficients.  The  bandwidth  for  each  region  will  be  — (for  the 


20kHz  and  30kHz  case  it  is  — )  and  this  can  be  a  limitation  for  realizing  a  filter  when  M 

6 

is  a  large  number  with  a  finite  number  of  sampling  points.  For  the  3  PRF  case,  30kHz, 
40kHz  and  50kHz,  the  total  number  of  filters  in  Figure  3.4.2  will  be  1 1  and  the  maximum 

i  i  11 

bandwidth  of  the  signal  would  be  \co\  <  k—  when  the  sampling  rate  is  600kHz  (=  l.c.m  of 

60 


30kHz,  40kHz  and  50kHz).  This  is  equivalent  to  sampling  the  data  at  110kHz  (= 
600x11/60)  since  there  is  both  a  reduction  of  the  bandwidth  and  an  increase  of  the 
sampling  frequency.  That  is,  the  maximum  resolvable  Doppler  frequency  is  increased  by 
2.75  (=110/40)  times  over  the  single  PRF  case  of  40kHz.  In  this  case  the  maximum 
resolvable  range  increases  by  4  (=40/10)  times  that  of  the  single  PRF  system  of  40kHz. 

2k  2k 

The  bandwidth  of  each  region  for  the  OMF  is  —  =  —  . 

M  60 


3.4.3  Optimum  Value  of  M  and  N  in  Radar  Application 

This  method  has  two  practical  limitations  as  indicated  earlier  which  causes  errors  in  the 
estimation.  First,  assume  that  the  signal  is  a  band  limited  sequence.  In  a  real  situation 
obtaining  a  band  limited  sequence  through  a  finite  number  of  samples  is  not  possible.  In 
radar  applications,  typically  the  maximum  number  of  samples  does  not  exceed  a  few 
hundred  points.  In  which  case,  our  reconstruction  of  spatial  domain  data  will  have  inherent 
errors.  This  is  shown  in  the  previous  example  of  Figure  3.4.5.  The  other  limitation  is  the 
bandwidth  of  each  region.  Each  of  the  filters  to  be  implemented  in  each  of  the  region  has 
an  M  band  filter  characteristic  [25]  and  the  bandwidth  becomes  too  small  if  M  becomes 
large.  For  the  3  PRF  case  of  30kHz,  40kHz  and  50kHz,  we  need  a  filter  which  has  a 
2k  2k 

bandwidth  of  —  =  —  .  This  type  of  filter  is  difficult  to  implement  since  there  are  only  a 
M  60 

few  hundred  sample  points. 
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The  range  increment  by  using  a  multiple  PRF  system  is  proportional  to  1/PRF  and  it  is 
given  by 


IR 


PRF2 

g.c.d {PRFl ,  PRF2  ■  ■  ’ 


(3-4-10) 


where  IR  is  the  increase  of  the  scale  factor  in  tenns  of  the  range  and  PRF1  is  the  PRF  used 
for  comparison  with  that  of  a  single  PRF  system. 

The  maximum  resolvable  range  for  the  30kHz,  40kHz  and  50kHz  PRF  cases  will  be  4 
(=  40/10)  times  that  of  the  single  PRF  case  corresponding  to  40kHz. 

The  increase  in  the  range  of  the  Doppler  by  using  a  multiple  PRF  system  can  be  seen 
from  the  following  equation 


ID  -  Bandwidth  x  j\ 


N  l.c.m(PRFl ,  PRF2  ■  ■  •) 
M  PRF2 


(3-4-11) 


where  ID  is  the  increase  in  the  range  of  the  Doppler,  Bandwidth  is  the  ratio  of  the 
maximum  signal  bandwidth  divided  by  two  times  the  Nyquist  rate  and  fs  is  the 
equivalent  sampling  frequency  of  the  multiple  PRF  system.  The  maximum  resolvable 
Doppler  for  the  30kHz,  40kHz  and  50kHz  PRF  case  is  1 1/4  times  that  of  the  one  PRF  case 
of  40kHz. 

Since  there  is  a  trade-off  between  maximum  Doppler  and  maximum  range  [1,  2],  we 
want  to  maximize  the  incremental  product  of  the  of  Doppler  and  the  range,  that  is, 


IDxIR 


N  LCM{PRFx,PRF2  •••)  PRFong 

M  ~PRF~  X  GCD(pRF1  ,  PRF2  •  •  •) 


_ n _ lcm{prfx,prf2 •••)  PRFolig 

lcm(prf1,prf2 •••)  PRFong  gcd{prfx,prf2 •••) 

gcd{prf„prf2 •••) 


=  N.  (3-4-12) 

We  can  maximize  N  with  respect  to  M.  From  this  result  one  can  choose  N  =  M- 1  since  N 
should  always  be  less  than  M.  For  the  case  of  M  =  16  and  Af  =15,  the  results  are  shown  in 
Figure  3.4.8  and  3.4.9  when  applied  to  the  same  signal  as  in  the  previous  example  (3-4-7). 
Note  that  Equation  (3-4-10)  is  true  for  the  general  multiple  PRF  system  and  Equation  (3- 
4-1 1)  is  applicable  only  when  using  the  QMF  filters.  Therefore  the  result  of  (3-4-12)  may 
be  different  if  other  reconstruction  techniques  are  used. 
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(a)  Spectrum  of  the  original  sequence 


(b)  Spectrum  of  the  estimated  sequence 


Frequency  [rad] 

Figure  3.4.8  Frequency  domain  response  for  the  optimized  case  (M  =  16,  N  =  15). 


(a)  Original  sequence 


Figure  3.4.9  Result  for  the  example  (magnitude)  in  the  time  domain  (M=16  and  N=\  5). 

3.5  Iterative  Method 

The  application  of  iterative  methods  based  on  the  Sandberg’s  theorem  [20]  was  first 
used  for  evaluating  the  spectrum  of  a  nonuniformly  sampled  data  by  Willy  [29].  Willy 
used  pulses  of  finite  width  and  conjectured  that  the  interactive  recovery  procedure  will 
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converge  for  the  finite  width  sampling.  Marvasti  [30]  tried  to  prove  that  the  iterative 
method  converges  for  random  samples  when  the  samples  are  chosen  from  a  unifonn  or  a 
Poisson  distribution.  Sandberg  proved  the  convergence  of  this  method  [31]. 

It  is  shown  that  if  the  signal  is  band-limited  and  the  average  sampling  rate  satisfies  the 
Nyquist  sampling  rate,  then  the  signal  can  be  recovered  [12,  13]  without  aliasing.  The 
method  uses  an  iterative  procedure  that  requires  low  pass  filtering  of  the  unequally  spaced 
samples  followed  by  resampling  at  the  same  points.  Again,  low  pass  filtering  is  applied  to 
the  signal  to  obtain  a  corrected  signal.  Repeated  application  of  this  process  is  shown  to 
converge  to  the  original  signal.  The  proof  of  convergence  is  given  in  Appendix  C.  By  its 
nature,  the  iterative  method  is  more  time  consuming  than  the  other  methods  since  it 
requires  a  large  number  of  iterations.  Even  though  the  convergence  is  guaranteed  by  this 
method,  the  number  of  iterations  can  be  very  large  to  reach  a  certain  error  criteria.  Some 
effort  has  been  done  by  Park  [24]  to  increase  the  rate  of  convergence  by  minimizing  the 
energy  term,  but  in  their  approach  one  needs  to  know  the  values  for  the  error  tenns  that  is 
usually  not  known.  Usually  a  smaller  spacing  between  samples  provides  a  faster 
convergence  of  the  sequence. 

Starting  with  the  Sandberg  theorem,  “On  the  properties  of  some  systems  that  distort 
signals“  [31]  Sandberg  extended  Beurling’s  theory  of  recovery  of  distorted  band-limited 
signals  in  a  Hilbert  space.  He  showed  that  signal  recovery  can  be  extended  to  the  case  in 
which  a  known  square  integrable  noise  is  added  to  the  input  signal  and  the  result  applied 
to  a  time  varying  device  which  may  be  non-linear. 

Sandberg’s  theorem  can  be  stated  as:  Let  P  and  Q  be  the  mapping  of  k  into  H  (Hilbert 
space)  such  that  for  all  f,ge  k 


Refer  -  Qg,  f  -  g)  2  K  \\f  -  g|f .  (3-5-1 ) 

\pQf-PGgf<k2\f-g\\  (3-5-2) 

where  kx  and  k2  are  positive  constants.  Then  for  each  h  e  k,  the  equation  h  =  PQf 
possesses  a  unique  solution 

( PQ )  '  h  e/c  ,  (3-5-3) 

given  by 

{PQ)  1  h  =  lim/n ,  (3-5-4) 

n— >oo 

where 

=  0-5-5) 

ft  2 

and  f0  is  an  arbitrary  element  of  k.  Furthennore,  for  all  h]  ,h2  e  k 


(pqY'K  -( PQYK 


VI 

h ,  -  h , 

1  k, 

1  z 

(3-5-6) 
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The  proof  is  given  in  Appendix  C.  The  uniqueness  of  the  solution  for  h  =  PQz  is  also 
given  in  Appendix  D. 

Applying  theorem  (3-5-6)  to  the  unequally  spaced  signal,  and  letting 
P  =  band-limiting  operator  (Low  pass  filter), 

Q  =  ideal  nonunifonn  sampling  operator  the  same  as  in  the  given  nonuniformly 
spaced  samples, 

will  transfonn  the  expressions  in  (3-5-4)  and  (3-5-6)  as 

\imxk{t)  =  x(t) ,  (3-5-7) 

k—>co 

xk+l(t)  =  XPQx(t)  +  (p-  XPQ)xk(t) ,  (3-5-8) 

where  x(t)  is  the  band-limited  signal  and  A,  is  a  constant.  The  constraint  between  the 
coefficient  X  and  data  spacing  t.  are  given  in  Appendix  F  and  is  given  by 

2k, 

0  <  2  <  — -,  (3-5-9) 

k2 

k2  <  k2  ■  (3-5-10) 

This  iterative  method  will  converge  provided  if 

k+1  (0-^(0  NlkW-^-lW||  for  all  k.  (3-5-11) 

We  have 


P{xk  {t)-xk-\  (t )) -  XPQ  {xk  (t )~xk_ |  (t ))  <r\xk{t)-xk_x{t)  ||  for  0  <  r  <1 .  (3-5-12) 


The  proof  of  (3-5-12)  is  given  in  Appendix  E.  Note  that  to  minimize  r  and  maximize  the 

k] 

rate  of  convergence,  it  is  required  that  X  =  — .  The  main  problem  in  choosing  a  proper  X 


is  that  it  cannot  be  detennined  theoretically  since  kx  and  k2  are  not  known.  X  is  generally 
chosen  between  0.5  and  1. 

Initial  values  of  the  analog  signal  x(t)  can  be  set  to  x(tk )  with  sample  and  hold,  which 


would  be  a  good  approximation  for  an  arbitrary  initial  value.  To  simulate  the  band-limited 
operator  Q  using  the  nonunifonnly  spaced  data,  the  FFT  has  been  used  using  the 
nonuniformly  spaced  data,  which  may  not  yield  an  exact  frequency  response  but  an 
approximation  to  it.  The  cut  off  frequency  with  zero  padding  has  been  set  to  0.6  n  ( n 
provides  an  all-pass  filter)  since  the  average  spacing  for  our  example  does  not  exceed  one 
half  of  the  Nyquist  sampling  interval.  Finally  the  inverse  FFT  is  used  to  get  the  band- 
limited  signal.  Here,  P  is  the  sampling  operator  which  converts  the  analog  signal  to  the 
nonuniformly  sampled  signal.  To  simulate  P,  the  Lagrange  polynomial  interpolation  has 
been  used  to  represent  the  function  and  is  evaluated  at  the  nonuniform  sampling  points. 
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One  constraint  of  the  iterative  method  is  the  sampling  points  tk  cannot  exceed  the 
range  of  the  A-th  interval  to  guarantee  convergence  [31],  that  is 


(3-5-13) 


where  Tk  is  the  location  of  the  sample  points  using  an  average  interval,  and  A  7)  is  the 

deviation  associated  with  the  A-th  interval  from  the  uniform  average  spacing. 

For  example,  consider  the  signal  of  (3-4-7)  which  is  given  by 

x{tk)  =  e2'2M*i  +2e-3-5-2^'  -2.5e4'2^';  for  k=  1,  2,  ...,  100, 


where  tk  is  the  time  step  satisfying  (3-5-13).  Let  the  deviations  from  the  uniform  spacing 
represent  a  set  of  uniformly  distributed  random  numbers  bounded  by  [-  ^ ^  ]  for  A  = 
1,  2,  ...,  100  and  T  =  0.03.  The  ratio  of  the  average  sampling  rate  to  the  Nyquist  sampling 
rate  is  33.33,/  Using  (3-5-8)  iteratively,  the  time  domain  sequence  and  the 

corresponding  spectrum  can  be  obtained  as  shown  in  Figure  3.5.1.  A  typical  plot  of  the 
iteration  number  as  a  function  of  the  error  is  given  in  Figure  3.5.2.  Theoretically  this 
method  has  been  shown  to  converge  but  the  error  does  not  go  to  zero  for  the  real  case. 
Here,  the  error  is  defined  as  the  normalized  difference  in  the  time  domain  data. 


(a)  Time  sequence 


Figure  3.5.1  Result  of  the  iterative  method 
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Figure  3.5.2  Reduction  of  error  with  iteration. 


3.6  Orthogonal  Polynomial  Expansions 

In  this  approach  the  unevenly  spaced  data  is  approximated  by  a  set  of  orthogonal 
polynomials.  Then  since  the  Fourier  transfonn  of  the  orthogonal  polynomials  are  known 
analytically  one  can  obtain  the  spectrum  of  the  unequally  spaced  data.  In  this  section  the 
Associated  Flennite  functions  and  the  Spherical  Bessel  functions  have  been  utilized  to 
recover  the  unevenly  spaced  signal  [32,  33].  The  Fourier  transforms  of  the  Associated 
Hermite  functions  and  the  Spherical  Bessel  functions  are  the  Associated  Hennite  functions 
and  the  Legendre  polynomials,  respectively. 

The  benefit  of  using  orthogonal  functions  is  that  they  are  easy  to  generate  by  using 
their  recursive  properties.  But  this  method  also  depends  on  the  spacing  between  samples. 
Moreover  there  is  no  measure  of  error  as  to  how  well  it  fits  the  original  signal.  Even 
though  the  fit  looks  good  in  the  spatial  domain,  there  is  no  guarantee  that  it  has  better 
estimation  properties  in  the  frequency  domain.  It  is  easily  seen  that  the  average  spacing 
should  be  much  less  than  the  Nyquist  sampling  rate  since  the  Nyquist  rate  allows  only  two 
samples  (or  more)  per  one  period  (for  the  case  of  a  real  signal).  Next,  the  formulation  in 
terms  of  orthogonal  polynomials  is  presented. 

3.6.1  Approximation  of  unevenly  spaced  data  by  the  Associate  Hermite  Polynomials 

The  associate  Hermite  functions  ( hn )  and  the  Hermite  functions  ( Hn )  are  related  by 
[14] 


K(t,i)  = 


V2"fl!  yj  -\l~7rl 


(3-6-1) 
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where  /  is  a  scaling  factor  and  n  is  the  degree  of  the  polynomial.  The  Hermite 
polynomials  can  be  computed  recursively  through 


H0{t)  =  1, 

Hl{t)  =  2t,  (3-6-2) 

Hn{t)  =  2tHn_x{t)-2(n-\)Hn_2(t)  ;  n>  2. 


Using  (3-6-1)  and  (3-6-2),  the  associated  Hennite  functions  can  be  calculated  easily 
through  the  recursive  relationship 


K(x) 


n>  2 . 


(3-6-3) 


Some  of  the  lower  degrees  of  the  associate  Hermite  polynomials  are  shown  Figure  3.6.1. 
Any  signal  can  be  expanded  by  the  associate  Hermite  functions  through 


where  hn 


™  n 
f{x,  )  =  Z  «  A  (*,-  )  =  Z  7 =K 

n= 0  n= 0  v  * 

4ihn[xi,i) . 


U 
V  i 


(3-6-4) 


Figure  3.6.1  A  plot  of  the  associate  Hennite  functions  of  different  degrees. 


Using  the  following  Equations  [14], 
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1 2  /  I  ^  00  _^,2  / 

e  r^H2m^)  =  (- !)'% -j e  ^  H2m(y)cosMdy’  (3-6-5-a) 

*  n  0 

^2/  (2”00  _^2  / 

e  t/2H2m+1(t)  =  (-l)mJ-fe  ^//2ra+1(v)sin(^y)fl?v,  (3-6-5-b) 

»  n  0 

a  Fourier  transform  relationship  by  combining  Equations  (3-6-1)  and  (3-6-5)  is  established 
as  [33] 


F(/)=ixw)" 


4hin 


-1 


2n+\ 


■\f^2 


J2»+l  \/l 


Y 1  t  -V'  /  f / 

=  ~n=hn\/ 

n= 0  2  / 


or  equivalently  the  Fourier  Transform  pair  is 


(3-6-6) 


^W(-0"4 


V7i 


(3-6-7) 


where  /2  =  /  is  a  scale  factor.  Thus,  if  we  can  expand  the  function  /(x;  j  by  an 

orthogonal  polynomial  hn ,  then  the  Fourier  transform  of  f{xj)  can  be  expressed  by 

with  the  same  coefficients.  The  coefficients 

2' 

an  can  be  calculated  from  the  matrix  representation 


adding  up  the  tenns  (-/)"  ~^=hn(^/i 


K{Xx) 

h2{x 1)  • 

••  ^(*1) 

a\ 

Ax  J 

K{xi) 

h2(x2)  ■ 

••  hN(x2 ) 

a2 

= 

Axi) 

K{xm) 

h2(xM)  • 

_aN  _ 

_Axm)_ 

where  M  is  the  number  of  data  points.  N  is  the  maximum  degree  of  the  associated 
Hermite  functions  and  hn  is  the  associated  Hermite  polynomial  of  degree  n.  Finally, 

Axj]  is  the  sampled  value  of  the  data  at  jcy  for  j  =  1,  2,  ...,  M.  The  degree  of  the 

Hermite  polynomial,  N,  can  be  set  up  to  a  maximum  value  which  is  identical  to  the 
number  of  data  samples  M.  Here,  x  -  is  not  limited  to  evenly  spaced  samples  only  but  also 

covers  unevenly  spaced  samples.  The  average  sampling  rate  needs  to  be  less  than  the 
maximum  signal  frequency  so  that  it  can  approximate  f(x)  in  a  smooth  fashion. 

As  an  example,  consider  the  signal  in  (3-4-7) 
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f(tk )  =  e2'27ttki  +  2e~3'5'2Mki  -  2.5e4'2a#*! ;  for  k  =  1,  2,  . . 100, 

where  tk  follows  a  uniformly  distributed  random  number  in  the  range  [-1,  1].  The 
average  sampling  rate  is  50  Hz  and  the  scale  factor,  /, ,  has  been  chosen  to  be  0.2.  an ’s 

have  been  obtained  using  (3-6-8)  and  the  corresponding  spectrum  can  be  estimated 
through  (3-6-6).  Figure  3.6.2-(a)  is  the  time  domain  signal  and  the  interpolated  function. 
Figure  3.6.2-(b)  is  the  corresponding  frequency  domain  response  using  the  associated 
Hermite  functions,  and  the  three  frequency  components  can  be  clearly  observed  in  the 
analyzed  data. 

3.6.2  Approximation  by  the  Legendre  Polynomial 

The  Legendre  polynomials  are  also  orthogonal  polynomials  and  their  Fourier 
transforms  are  analytically  known,  which  are  the  spherical  Bessel  functions.  The  Legendre 
polynomials  exist  only  in  the  region  -1  to  1.  A  complex  function  of  finite  support  can  be 
expanded  in  terms  of  these  polynomials,  as  the  time  domain  sequence  is  finite.  The  band- 
limited  Legendre  functions  can  be  scaled  to  any  arbitrary  range  for  practical  use.  The 
relationship  between  the  spherical  Bessel  functions  and  the  orthogonal  Legendre 
polynomials  is  given  [14]: 

P x~^e~imJ  u  {x)dx  =  (-;)"( 2jc)& pAcq)  ;  co2  <  1  (3-6-9) 

J— oo  w+/2 

=  0  ;  CO2  >  1. 

(a)  Time  sequence 
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(b)  Estimated  spectrum 
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Figure  3.6.2  Fitting  of  the  data  by  the  associate  Hermite  functions. 
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Spherical  Bessel  functions  can  be  written  in  terms  of  the  ordinary  Bessel  functions  as: 


Jn(z) 


(3-6-10) 


Substituting  (3-6-10)  into  (3-6-9)  will  yield 

|  e  im  — -  j  (x)dx  =  Pn(co)  ;  \co\  <  1 .  (3-6-11) 

J-CO  n 

Using  the  Fourier  transform  property, 


l*co  ^ 

p  (x)  =  j  2(-  i)  j  (-  co)e  'mdco  ;  - 1  <  x  <  1 

»  — CO 


(3-6-12) 


and  by  using  a  scale  factor  a,  one  obtains 

Pn^— j  =  1 2{- i)n\a\j n{- aco )e~,aK da>  ;  -  a  <  x  <  a  .  (3-6-13) 

^  —co 

n 

Different  orders  of  Pn(x )  and  2 (-/')  jn(  -  co)  are  shown  in  Figures  3.6.3(a)  and  (b).  For 
example,  the  rectangular  function  which  corresponds  to  a  Legendre  function  with  n= 0 
whose  Fourier  transform  corresponds  to  a  spherical  Bessel  function  of  zero  degree  which 
is  the  sine  function  as  shown  in  Figure  3.6.3(b).  Therefore,  if  we  can  expand  the  function 
f(x)  by  using  the  orthogonal  Legendre  polynomials  Pn ,  then  the  Fourier  transform  of f(x) 

can  be  expressed  by  adding  up  the  spherical  Bessel  functions  with  the  same  orders.  Or,  if  a 
signal  can  be  expanded  by  the  spherical  Bessel  functions  then  the  Fourier  transform  can  be 
evaluated  using  the  corresponding  sum  of  Legendre  polynomials. 

A  signal,  for  example,  can  be  written  in  terms  of  an  infinite  sum  of  Legendre  functions 
as 


f{X)  =  HAnPXX)' 


n= 0 


(3-6-14) 


Multiplying  both  sides  by  Pm\Xj)  and  integrating  from  -  oo  to  oo  will  give 


OO  OO  00 

I  f(xi )  Pm  (x{  )dx  —  |  X  dn  Pn  (xj )  Pr)l  (xj )  dx 

—oo  —  oo/I^O 

oo  oo 

=  I  1  An  Pn(xi)Pm(xi)  dx  ■  (3-6-15) 

n= 0  —oo 
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Real  part  Imaginary  part 


(b) 

n 

Figure  3.6.3  (a)  Legendre  functions  of  different  degrees  and  (b)  The  plot  of  2(—  i)  j n  (—  to  )  for  different 

degrees 


The  right  hand  side  will  be  zero  except  for  m  =  n  since  P j ’s  are  orthogonal  to  each 
other  and  the  left  hand  side  of  the  integral  will  then  be  a  summation. 

Z/hhh)=S4.[^h)]2  ■  (3-6-16) 


Therefore  the  coefficients  An  are  given  by 
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(3-6-17) 


An  =Z 


f{Xj)P»{Xj) 


~  r„ 

N  2 

where  yn  =  X  (  ^  (x/  ))  •  N  's  the  number  of  data  and  Pn 

j= i V  ' 

degree  n.  Then, 


f(^)  =  TiA„  2(-  1)"a(®). 


is  the  Legendre  polynomial  of 


(3-6-18) 


Note  that  the  Bessel  functions  are  defined  only  along  the  positive  axis. 

Since  the  transform  is  limited  from  -1  to  1  (or  -a  to  a  through  scaling),  the  frequency 
response  would  consist  of  sine  functions  which  will  produce  undesired  aliasing  between 
the  various  frequency  components.  To  obtain  a  better  estimate  of  the  amplitude,  the  same 
technique  as  in  the  previous  section  can  be  used. 

As  an  example,  consider  a  signal  as  defined  by  (3-4-7) 

f(t  k)  =  e22Mti  +  2e-3-5'2^'  -2.5e4'2^';  for  k=  1,  2,  ...,  100, 

where  tk  is  generated  for  a  uniformly  distributed  random  number  in  the  range  of  [-1;  1]. 
The  average  sampling  rate  is  50Hz.  An’s  have  been  obtained  using  (3-6-17)  and  the 
corresponding  spectrum  can  be  estimated  from  (3-6-18).  Figure  3.6.4(a)  is  a  plot  of  the 
time  domain  signal  and  the  reconstructed  signal.  Figure  3.6.4(b)  is  the  corresponding 
frequency  domain  response  using  the  spherical  Bessel  functions.  Note  that  the  spherical 
Bessel  functions  only  exist  for  the  positive  argument  and  thus  the  negative  frequency  in 
the  signal,  3.5  Hz,  cannot  be  represented  in  the  Figure  3.6.4(b). 


3.7  Estimation  in  terms  of  the  Analog  Frequency 

One  of  the  techniques  of  performing  a  perfect  reconstruction  of  a  band-limited  signal 
from  its  unevenly  sampled  data  points  has  been  proposed  by  Jenq  [34,  35,  36].  This 
method  shows  that  if  the  average  sampling  rate  is  greater  than  the  half  of  the  maximum 
frequency  component  of  the  signal,  then  the  signal  can  be  perfectly  reconstructed.  While 
the  Lomb  periodogram  is  not  affected  very  much  by  the  sampling  rate  since  it  estimates 
the  amplitude  and  phase  by  a  Least  squares  technique,  Jenq’s  and  the  QMF  methods 
produce  perfect  reconstruction  from  a  theoretical  point  of  view  if  the  average  sampling 
rate  satisfies  the  Nyquist  rate.  In  comparison  to  the  QMF  method,  Jenq’s  method  does  not 
need  to  use  a  set  of  filter  banks  which  can  cause  some  errors.  Here  the  spacing  between 
the  data  samples  cannot  be  random.  That  is  the  deviation  from  a  uniformly  sampled  data, 
denoted  by  rm ,  should  be  limited  in  magnitude  as  in  the  case  of  the  QMF  method.  The 
nonunifonn  sampling  rate  can  be  expressed  as 

t„=nT  +  A„,  (3-7-1) 
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where  T  is  the  period  of  the  uniformly  sampled  signal.  A((  is  the  periodic  sequence  with  a 
period  M.  Let  n  =  kM  +  m,  where  k  ranges  from  -oo  to  +oo  and  m  ranges  from  0  to 
(M-l),  then 


tn  =  {kM  +  m)T  +  A/h 

=  kMT  +  mT  +  rmT  ,  (3-7-2) 


where  rm  =  ^m/j  is  the  timing  offset  expressed  in  percentage  in  terms  of  the  nominal 
sampling  period  T.  The  digital  spectrum  for  the  nonuniform  samples  is  given  by 


oo 

xd(<°)=  Zx(0  e~JU"  ■ 

k=— oo 


(3-7-3) 
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Figure  3.6.4  Example  using  the  Legendre  polynomial 


Utilizing  (3-7-2)  will  yield 


oo  M- 1 

Xd{(0)=  Y,  YX(kMT  +  mT  +  rmT)  e-MkMT+mT+r-T)  . 

k=— oo  m= 0 


(3-7-4) 
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Since  x(kMT  +  mT  +  rj)  =  —  [  Xa  (q)  e/n(fcwr+mr+r"'r)JQ ,  one  obtains 

2  n J~c0 


oo  M-l  1 

*»  =  Z  e'n("r*"r"-rl<<n  e 

M=n 


j£l{kMT +mT +rmT)  —jcot(kMT +mT +rmT ) 


k=-cc  m= 0  - 


M-l  1  00 

-li  *.(n) 


z< 


J(Cl-co)kMT 


k=-  co 


j{Cl-o)){mT +rmT ) 


dCl , 


(3-7-5) 


m=0  2^r 

and  the  following  discrete  sequence  can  be  obtained  by  using  the  sum  of  delta  functions  as 


M-l  i  °° 

m=o2x  L 


2n 


r 


Y—8\ 

kTLMT  1 


Q-<z>  +  £ 


V 


2;r 

MT 


j  (Cl-co)(mT +rmT ) 


r/Q . 


(3-7-6) 


By  changing  the  order  of  integration  and  summation,  one  obtains 


oo  M-l 


i=-t»  m=0  lvl  ‘  _ qq  V 


Q  -  co  +  k 


2n 

MT 


j  (Q-(o)(mT +rmT ) 


dCl 


1  rtxL-^ 


MT 


k=- oo  m=0 


MT 


3  -jk - [mT+r  T) 

e  MT 


1  V-I  1  X— I1  -jk'm—  -jkm— 

-  Y  —Ye  Me  M 

Tt\Mt« 


In  Tn  A 


\ 

(  2  n  ^ 

co-k - 

) 

v  mt) 

(3-7-7) 


Finally, 


1  7 

^(®)=-Zuo^ 

1  k=- oo 


f  2n^ 

co-k - 

l  MT) 


(3-7-8) 


where 


1  -/fo3 

Y(k)  =  —  ye  Me  M 

M  m=0 


(3-7-9) 


(3-7-8)  shows  the  relationship  between  the  discrete  spectrum  and  the  analog  spectrum.  If 
we  consider  the  angular  frequency  between  zero  to  2  n  then 


*A<»)  =  k  tr(k)x_ 


T  M 

k= - +1 

2 


co-k 


2  n 
MT 


(3-7-10) 


To  find  the  analog  spectrum  X a ,  one  can  change  the  variable  k  from  0  to  M-l.  That  is, 
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X.. 


co  +  i  n  - 


V 


2  n 
MT 


7  2  Ak  +  m)  X, 


T  M 

k= - +1 

2 


co-k 


2n 

MT 


(3-7-11) 


or  equivalently 


[x»]  =  [r]  [*»], 


(3-7-12) 


where 
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X. 
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xdM 
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MT 
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(3-7-13) 


By  taking  the  inverse  of  the  M  by  M  matrix  Y  and  multiplying  it  by  Xd  one  can  get  the 
analog  frequencies  Xa . 

To  illustrate  the  applicability  for  the  multiple  PRF  system,  we  take  the  signal  of  (3-1- 
2),  i.e., 


f{xk)  =  e22mki  +2e-3S2mti  - 2.5e4'lMXki ;  for  k=  1,2,...,  128. 

The  signal  is  sampled  at  two  different  frequencies  simultaneously.  The  two  sampling 
frequencies  are  PRFX  =  6.25Hz,  PRF2=9.375Hz.  This  is  equivalent  to  sampling  by  2Hz 
and  3Hz  except  for  the  scaling  in  frequency.  PRFX  has  2  samples  in  a  record  of 
0.32seconds,  PRF2  has  3  samples  in  0.32seconds  and  this  repeats  every  0.32  seconds. 
This  will  occur  every  0.32seconds  and  the  total  number  of  samples,  M,  is  obtained  by 
adding  up  the  samples  for  all  the  PRFs.  Since  the  last  sample  in  the  record  of  0.32 
seconds  is  also  the  starting  sample  for  the  next  time  interval,  M  will  be  given  by  M  — 
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PRFX  +  PRF2  -  1,  if  the  PRFs  are  relative  prime  numbers.  The  matrix  Y  can  be  computed 
from  (3-7-9)  and  (3-7-13)  and  is  used  to  obtain  Xd(co).  The  result  of  this  approach  is 
compared  to  that  of  the  DFT  of  the  evenly  spaced  data  as  shown  in  Figure  3.7.1  and  3.7.2. 
Figure  3.7.3  is  the  DFT  of  the  unevenly  spaced  signal.  It  is  seen  that  the  weak  signal 
cannot  be  distinguished  in  Figure  3.7.3  but  can  be  discerned  quite  easily  in  Figure  3.7.2. 


Frequency  [Hz] 


Figure  3.7.1  DFT  of  the  uniformly  spaced  signal. 


Frequency  [Hz] 


Figure  3.7.2  Spectrum  of  the  nonuniformly  sampled  signal  after  processing 
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25 


Frequency  [Hz] 

Figure  3.7.3  DFT  of  the  nonuniformly  spaced  data. 


3.8  Comparison  of  the  Various  Methods 


The  performance  of  all  the  methods  described  so  far  is  now  compared  with  each  other 
in  this  section.  Since  for  every  method  presented,  there  is  a  problem  of  some  sort 
associated  with  each  technique,  it  is  difficult  to  compare  the  results  from  all  the 
approaches  directly.  Usually,  the  polynomial  type  methods  as  described  in  section  3.1  and 
section  3.6  do  not  depend  much  on  the  number  of  data  samples.  Rather  they  fit  the  data 
well  with  a  less  number  of  points  because  the  degree  of  the  polynomial  becomes  large  and 
the  error  associated  with  the  approximation  of  the  data  may  become  great  due  to 
numerical  instabilities.  The  other  methods  perform  better  as  the  number  of  data  samples 
increases.  The  QMF  method  (section  3.4)  needs  more  than  a  few  hundred  data  points  to 
result  in  a  meaningful  estimation.  The  use  of  orthogonal  polynomials  and  the  iterative 
method  need  to  have  a  small  sampling  interval  while  the  Least  squares  approach  (of 
section  3.3)  is  not  much  affected  by  the  average  sampling  rate.  We  divide  all  the  examples 
into  three  parts.  We  consider  the  following  three  cases,  (1)  in  this  case,  the  sampling 
frequency  is  much  greater  than  the  highest  frequency  of  the  signal,  (2)  the  sampling 
frequency  is  equal  to  the  highest  frequency  of  the  signal  and  finally,  (3)  the  sampling 
frequency  is  much  lower  than  the  highest  frequency  of  the  signal  in  which  case  aliasing 
may  occur.  Note  that  the  sampling  frequency  defined  in  this  context  to  recover  the  original 
complex  signal  is  equal  to  the  maximum  frequency  of  the  signal.  It  is  not  related  to  the 
conventional  Nyquist  sampling  frequency.  For  each  case,  the  results  computed  by  the 
various  methods  are  compared.  To  illustrate  the  accuracy  and  the  stability  of  the 
techniques,  ground  clutter  has  been  added  in  addition  to  white  Gaussian  noise  and  the  final 
result  is  compared  with  that  from  the  conventional  CR  theorem.  For  all  the  examples,  we 
considered  a  2  PRF  system  using  7Hz  and  8Hz  as  the  two  sampling  frequencies.  The 
average  sampling  frequency  is  14Hz  since  it  has  14  samples  within  one  time  step. 
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According  to  (2-10),  the  maximum  resolvable  Doppler  frequency  is  increased  to  8  times 
as  compared  to  the  case  of  a  single  PRF  system  of  7Hz.  The  number  of  sampling  points 
used  by  each  method  is  also  presented. 

3.8.1  Case  1:  fs  >  fmax 

As  mentioned  in  the  previous  sections,  the  polynomial  interpolation  method,  iterative 
method  and  orthogonal  expansion  method  yield  better  estimation  for  the  spectrum  when 
the  sampling  frequency  becomes  much  higher  than  the  Nyquist  frequency.  These 
approaches  usually  provide  a  meaningful  result  when  the  average  sampling  rate  is  at  least 
4  times  that  of  the  maximum  frequency  of  the  signal.  In  the  numerical  example,  the 
sampling  frequency  has  been  chosen  to  be  4  times  that  of  the  maximum  signal  frequency. 
For  the  first  example,  we  examine  a  single  frequency  signal  of  magnitude  10. 

Example  1:  Single  frequency  f  =  4/max 

The  average  sampling  rate  is  14Hz  and  the  signal  frequency  is  3.5Hz.  The  non-uniformly 
spaced  samples  repeat  every  14th  sample  and  a  total  of  100  sampling  points  are  taken  so 
that 


f{xk)  =  \0e2m'35Xk  with  k=\,2,...,N.  N  =  number  of  samples  =  100. 

The  nonunifonnly  sampled  time  domain  signal  and  the  DFT  of  the  original  signal  are 
shown  in  Figure  3.8.1(a)  and  (b).  The  frequency  response  has  a  peak  at  the  signal 
frequency  component,  which  is  the  25th  cell  out  of  100  cells.  All  the  eight  methods  have 
been  simulated  and  the  results  are  given  in  Table  3.8.1.  The  error  in  the  time  domain  has 
been  computed  by  using  the  normalized  average  of  the  differences.  The  estimate  of  the 
signal  is  f(xk )  and  the  error  in  the  approximation  of  the  function  is  defined  by 


Err  or  l  = 


N 


f(xk)~Hxk) 


N 


max 


\f(xk)\ 


(3-8-1) 


The  estimate  of  the  error  in  the  frequency  domain  is  given  in  Table  3.8.1.  It  is  computed 
by  taking  the  differences  between  the  maximum  of  the  computed  frequency  response  and 
the  signal  frequency  indicated  in  the  last  column  of  Table  8.3.1  and  is  expressed  by 


Error 2  = 


/-/ 


(3-8-2) 


where  /  is  the  target  frequency  and  /  is  the  estimated  frequency.  /  is  the  frequency  at 
which  the  maximum  value  of  the  frequency  response  is  obtained. 

Note  that  the  frequency  domain  error  tenn  is  the  difference  between  the  computed 
spectrum  of  the  time  domain  data  and  the  actual  frequency  response.  For  the  method 
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based  on  the  CR  theorem,  we  took  53  and  48  sample  points  for  each  sampling  frequency. 
Observe  that  the  peaks  in  the  spectrum  occurred  at  the  23ld  and  24th  cell  since  the  sampling 
frequency  of  8Hz  and  7Hz  will  yield  the  cell  number  23  (»  3.5 -53/8)  and  24 
(«  3.5  -48/ 7)  for  a  target  frequency  of  3.5.  Even  though  the  time  domain  estimation  is 
not  very  good,  acceptable  estimates  for  the  target  Doppler  can  be  obtained  for  all  the 
methods.  In  this  example,  all  the  methods  provided  acceptable  estimates  for  the  spectrum. 


(a)  Time  domain  sgnal 


Table  3.8.1  Summary  of  the  results  for  Example  1.  fs  =  4  fmax  with  single  frequency. 


Methods 

Errorl 

(time  domain) 

Error2 

(frequency  domain) 

Number 
of  data 

Target  cell 
number 

Lagrange  polynomial 

0.3077 

0.0000 

44 

11 

Cauchy’s  method 

0.9267 

0.0000 

44 

11 

CR  theorem 

- 

0.0000 

48,53 

24,23 

Least  squares  method 

- 

0.0000 

100 

25 

QMF  method 

0.2135 

0.0000 

100 

25 

Iterative  method 

0.3294 

0.0000 

100 

25 

Orthogonal  polynomial  approach 

2.6992E-4 

0.0000 

100 

25 

Analog  frequency  approach 

- 

0.0000 

100 

25 

Example  2:  A  Single  frequency  with  noise 


The  same  signal  as  in  Example  1  has  been  considered  and  additive  white  Gaussian  noise 
(AWGN)  has  been  added.  The  model  for  the  noise  signal  is  described  by 
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Noise  = 


(3-8-3) 


where  Pnoise  is  the  noise  power  and  N(0, 1 )  follows  a  normal  distribution  with  unit  variance 
and  zero  mean.  The  signal  model  is  defined  by 

f(xk )  =  lOe2®'3'5**  +  Noise  . 

Table  3.8.2  presents  the  probability  of  detection  for  each  method  as  a  function  of  the 
power  associated  with  the  noise  process.  We  assume  that  a  target  has  been  detected  when 
the  computed  spectrum  has  a  peak  which  corresponds  to  the  appropriate  cell  number  for 
the  target.  Total  number  of  trials  for  each  method  used  to  compute  the  probability  of 
detection  is  100.  With  an  increase  in  SNR,  the  probability  of  detection  increases  as 
expected. 


Table  3.8.2  Result  of  Example  2. 


fs  =  4  fm 


with  a  single  frequency  and  additive  AWGN 


Methods 

Probability  of  detection 

Number 
of  data 

Target  cell  number 

SNR  [dB] 

0  5  10  15  20  25 

Lagrange  polynomial 

0.97  1.00  1.00  1.00  1.00  1.00 

44 

11 

Cauchy’s  method 

0.18  0.27  0.58  0.70  0.99  1.00 

44 

11 

CR  theorem 

1.00  1.00  1.00  1.00  1.00  1.00 

48,  53 

24,  23 

Least  squares  method 

1.00  1.00  1.00  1.00  1.00  1.00 

100 

25 

QMF  method 

1.00  1.00  1.00  1.00  1.00  1.00 

100 

25  ! 

Iterative  method 

1.00  1.00  1.00  1.00  1.00  1.00 

100 

25 

Orthogonal  polynomial  approach 

0.00  0.00  0.00  0.00  0.00  0.00 

100 

25 

Analog  frequency  approach 

1.00  1.00  1.00  1.00  1.00  1.00 

100 

25 

Example  3:  Multiple  frequency  signals  and  f  =  4/max 
In  this  example,  the  signal  has  5  frequency  components 

^  =  \Oe2™035 Xk  +  ge2*105v  +  ge2,ri'lAxt  +  2elnil'iXk  +  6e2”3'5;tt 

We  assumed  that  the  last  component,  6e2m'2Sxt ,  is  the  signal  of  interest  (SOI)  and  the  other 
components  are  interference.  The  frequency  scanning  is  performed  from  2Hz  to  15Hz 
(equivalent  from  14th  to  100th  cell  in  the  frequency  domain  response)  since  the  major 
interfering  signals  are  located  below  2Hz.  The  time  domain  signal  and  the  corresponding 
frequency  domain  response  are  shown  in  Figure  3.8.2.  Table  3.8.3  is  the  result  of  the 
simulations.  The  CR  theorem  approach  has  been  deleted  from  the  Table  since  it  can  only 
deal  with  a  signal  having  a  single  frequency  component.  The  highest  degree  of  polynomial 
for  the  Cauchy’s  method  and  for  the  Lagrange  polynomial  technique  has  been  a  bit  higher 
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than  in  Example  1  since  for  Example  2  the  signal  has  a  stronger  low  frequency  component 
than  in  Example  1 . 

Example  4:  Multiple  frequency  signal  with  noise 

For  this  example,  the  signal  has  five  frequency  components  and  is  given  by 

■J=10e2*'°-35lt  +%e2M'lMXk  +  9e27lilAxt  +  2e27ri2Sxt  +  6e2”'3'5Xi  +  noise 

The  noise  model  is  the  same  as  in  (3-8-3)  and  the  SOI  is  the  same  as  in  Example  3.  The 
results  are  presented  in  Table  3.8.4  using  100  trials.  Except  for  the  Cauchy’s  method,  all 
the  other  methods  provided  acceptable  results  for  different  values  of  the  SNR.  Because 
Cauchy’s  method  provides  a  fit  to  the  data  through  the  use  of  rational  polynomials,  it 
cannot  precisely  fit  the  noisy  data  that  changes  rapidly  with  time. 


(a)  Time  domain  sgnal 


Figure  3.8.2  Sample  signal  for  Example  3 


Table  3.8.3  Summary  of  the  results  for  Example  3.  f  =  4  _/"  with  5  signal  components 


Methods 

Errorl 

(time  domain) 

Error2 

(frequency  domain) 

Number  of 
data 

Target  cell 
number 

Lagrange  polynomial 

0.0454 

0.0000 

44 

11 

Cauchy’s  method 

0.0755 

0.0000 

44 

11 

Least  squares  method 

- 

0.0000 

100 

25 

QMF  method 

0.0785 

0.0000 

100 

25 

Iterative  method 

0.0626 

0.0000 

100 

25 

Orthogonal  polynomial  approach 

7.1 196E-5 

0.0000 

100 

25 

Analog  frequency  approach 

- 

0.0000 

100 

25 
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Table  3.8.4  Summary  of  results  for  Example  4.  fs—  4  fmm  with  5  signal  components  and  AWGN 


Methods 

Probability  of  detection 

Number  of 
data 

Target  cell 
number 

SNR  [dB] 

0  5  10  15  20  25 

Lagrange  polynomial 

0.97  1.00  1.00  1.00  1.00  1.00 

44 

11 

Cauchy’s  method 

0.62  0.75  0.83  0.90  0.80  0.88 

44 

11 

Least  squares  method 

1.00  1.00  1.00  1.00  1.00  1.00 

100 

25 

QMF  method 

1.00  1.00  1.00  1.00  1.00  1.00 

100 

25 

Iterative  method 

1.00  1.00  1.00  1.00  1.00  1.00 

100 

25 

Orthogonal  polynomial  approach 

0.00  0.00  0.00  0.00  0.00  0.00 

100 

25 

Analog  frequency  approach 

1.00  1.00  1.00  1.00  1.00  1.00 

100 

25 

In  summary,  when  the  sampling  frequency  is  high  compared  to  the  highest  frequency 
content  of  the  signal,  most  of  the  methods  presented  in  this  report  can  be  applied  to 
estimate  the  spectrum  of  a  nonuniformly  sampled  data  sequence.  But  this  case  is  of  little 
use  to  the  radar  community  since  the  maximum  resolvable  Doppler  frequency  is  not  very 
large  as  seen  from  (2-10).  In  the  next  case  we  will  increase  the  signal  frequencies  (or 
decrease  the  average  sampling  rate)  and  will  observe  which  approach  has  a  better 
performance. 

3.8.2  Case  2:  fs  =  2/max 

When  the  sampling  frequency  equals  twice  the  value  of  the  maximum  signal 
frequency,  all  polynomial  methods  approximating  the  time  domain  samples  cannot  be 
used  since  they  need  smaller  sampling  periods.  Moreover,  the  results  of  the  QMF  method 
is  only  valid  for  the  case  when  the  average  sampling  rate  is  higher  than  the  Nyquist 
sampling  rate  and  its  results  will  degrade  in  accuracy  if  it  is  less  than  or  equal  to  the 
Nyquist  sampling  rate.  Only  three  methods,  the  Least  squares  method,  the  CR  theorem 
method  and  the  analog  frequency  approach  can  be  used  for  this  particular  case.  In  this 
example,  the  other  conditions  are  the  same  as  in  case  1  except  for  the  value  of  the  signal 
frequency. 


Example  5: 

The  signal  used  in  this  case  is  given  by  f{xk)  =  \0e2mlXk  and  the  results  are  given  in 

Table  3.8.5.  All  the  three  methods  yield  good  results.  (3-8-1)  and  (3-8-2)  has  been  used 
to  get  the  error  estimate. 


Table  3.8.5  Summary  of  the  results  for  Example  5.  fs  =  2  fnm  with  a  single  frequency 


Methods 

Error  1 

(time  domain) 

Error2 

(frequency  domain) 

Number  of 
data 

Target  cell 
number 

CR  theorem 

- 

0.0000 

48,  53 

0,46 

Least  squares  method 

- 

0.0000 

100 

50 

Analog  frequency  approach 

- 

0.0000 

100 

50 
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Example  6:  In  this  case,  the  signal  is  given  by  f{xk )  =  1 0e2"'’7x‘  +  Noise  and  the  noise  is 

computed  using  (3-8-3).  The  results  are  given  in  Table  3.8.6.  All  three  methods  yield 
good  results. 


Table  3.8.6  Summary  of  the  results  for  Example  6.  f  =  2  fnwx  with  a  single  frequency  and  AWGN 


Methods 

Probability  of  detection 

Number  of  data 

Target  cell 
number 

SNR  [dB] 

0  5  10  15  20  25 

CR  theorem 

1.00  1.00  1.00  1.00  1.00  1.00 

48,  53 

0,  46 

Least  squares  method 

1.00  1.00  1.00  1.00  1.00  1.00 

100 

50 

Analog  frequency  approach 

1.00  1.00  1.00  1.00  1.00  1.00 

100 

50 

Example  7:  In  this  case  the  signal  is  of  the  fonn 

f(xk)=l  0e2'n'OJXk  +8e2m'21xt  +  ge2*2-***  +  2e2^5-6^  +  6el7CilXk 

and  the  results  are  given  in  Table  3.8.7.  The  clustering  algorithm  was  not  used  in  this 
example  when  using  the  method  based  on  the  CR  theorem  since  there  are  many  frequency 
components,  and  therefore  only  the  Least  squares  and  the  analog  frequency  approach  can 
provide  acceptable  results. 


Table  3.8.7  Summary  of  results  for  Example  7.  fs  =  2  fmm  with  five  signal  frequencies 


Methods 

Error  1 

(time  domain) 

Error2 

(frequency  domain) 

Number  of 
data 

Target  cell 
number 

Least  squares  method 

- 

0.0000 

100 

50 

Analog  frequency  approach 

- 

0.0000 

100 

50  | 

Example  8:  In  this  example,  the  signal  plus  noise  is  described  by 

f(xk)  =  lOe2*'0'7'*  +  8e2™'21xt  +  9e2*i'Z8Xk  +2e2’ri5(,Xt  +  6e2”'7x 


‘  +  Noise 


and  the  results  are  given  in  Table  3.8.8.  The  clustering  algorithm  was  not  used  for  the 
method  based  on  the  CR  theorem  as  in  this  example  there  are  many  signals  at  different 
frequencies,  and  therefore  only  the  Least  squares  method  and  the  analog  frequency 
approach  yield  good  results. 

Example  9:  Multiple  signal  frequencies  with  added  noise  and  ground  clutter 

In  this  case,  the  signal  has  5  frequency  components  along  with  additive  Gaussian  noise 
and  ground  clutter. 

f(xk)=  I0e2”°-35XA  +%e2*isaxk  +9e2m'Axk  +2e2M2*Xk  +  6e2”'3'5^  +  noise  +  clutter. 
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Table  3.8.8  Summary  of  results  for  Example  8.  fs=  2  fmax  with  5  signal  frequencies  and  AWGN 


Methods 

Probability  of  detection 

Number  of  data 

Target  cell 
number 

SNR  [dB] 

0  5  10  15  20  25 

Least  squares  method 

0.97  1.00  1.00  1.00  1.00  1.00 

100 

50 

Analog  frequency  approach 

1.00  1.00  1.00  1.00  1.00  1.00 

100 

50 

The  clutter  model  is  given  by 

C(a)  =  2000  cT(01<B)2  (3-8-4) 

and  the  corresponding  frequency  spectrum  is  shown  in  Figure  3.8.3.  The  noise  model  is 
the  same  as  in  (3-8-3)  and  the  SOI  is  the  same  as  in  Example  3.  The  results  are  given  in 
Table  3.8.5  by  using  100  trials.  The  frequency  is  scanned  from  2FIz  to  7.5Hz  (equivalent 
from  14th  to  50th  cell  in  the  frequency  domain  response)  since  the  major  portion  of  the 
interference  and  the  ground  clutter  is  located  below  2Hz.  The  simulation  result  is  given  in 
Table  3.8.9. 

2000 
1500 
1000 
500 


0  10  20  30  40  50  60  70  80 

Frequency 

(b)  Received  signal  model  in  the  frequency  domain 

2000 
1500 
1000 
500 


0  10  20  30  40  50  60  70  80 

Frequency 

Figure  3.8.3  Model  for  a  ground  clutter  (evenly  sampled  data) 


Table  3.8.9  Summary  of  the  results  for  Example  9.  fs  =  2  fmax  with  5  signal  frequencies  along  with 

AWGN  and  ground  clutter. 


Methods 

Probability  of  detection 

Number  of  data 

Target  cell 
number 

SNR  [dB] 

0  5  10  15  20  25 

Least  squares  method 

1.00  1.00  1.00  1.00  1.00  1.00 

100 

50 

Analog  frequency  approach 

1.00  1.00  1.00  1.00  1.00  1.00 

100 

50 

(a)  Ground  clutter  model  in  the  frequency  domain 
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In  summary,  when  the  sampling  frequency  is  the  same  as  the  Nyquist  sampling  rate  of  a 
real  signal,  only  three  methods  can  be  applied  to  evaluate  the  spectrum  of  a  nonuniformly 
sampled  data.  In  this  case,  it  is  possible  to  achieve  a  higher  resolution  as  illustrated  in  (3- 
4-12).  By  using  the  three  methods  illustrated  in  this  case,  10-20  times  performance 
enhancement  in  terms  of  the  maximum  resolvable  Doppler  can  easily  be  obtained.  As 
mentioned  earlier,  the  method  based  on  the  CR  theorem  cannot  be  used  if  multiple  signals 
exist  or  if  there  is  ground  clutter.  In  the  next  example,  we  will  further  increase  the  number 
of  signal  frequencies  and  will  observe  which  approach  has  a  better  performance. 


3.8.3  Case  3:  fs  <  fmax 

When  the  average  sampling  frequency  is  less  than  the  maximum  signal  frequency,  only 
the  Least  squares  method  and  the  clustering  algorithm  can  be  used.  All  other  conditions 
remain  the  same  as  in  case  1  and  2  except  for  the  signal  frequency. 


Example  10:  For  this  example  fs  =  ^-/max  with  a  single  signal  frequency. 
The  data  is  given  by 


/(x.LlOe2"2 


and  the  results  are  given  in  Table  3.8.10.  All  the  three  methods  yield  good  results.  Here, 
(3-8-1)  and  (3-8-2)  have  been  used  to  obtain  an  estimate  of  the  error  in  the  time  and 
frequency  domain.  Note  that  the  error  in  the  clustering  algorithm  is  due  to  the  FFT  of  the 
data  samples,  which  do  not  match  exactly  to  the  signal  frequency. 


Table  3.8.10  Summary  of  results  for  Example  10.  fs  =  —  fmx  with  a  single  frequency 


Methods 

Error  1 

(time  domain) 

Error2 

(frequency  domain) 

Number  of 
data 

Target  cell 
number 

CR  theorem 

- 

4.4881E-4 

48,  53 

0,13 

Least  squares  method 

- 

0.0000 

too 

300 

Example  11:  In  this  case,  f{xk)  =  l0e2m'42Xk  +  Noise  and  the  noise  are  given  by  (3-8-3). 
The  results  are  given  in  Table  3.8.11.  Both  methods  yield  good  results. 


Table  3.8.11  Summary  of  the  results  for  Example  11.  fs  =  —  /max  with  a  single  frequency  and 

AWGN 


Methods 

Probability  of  detection 

Number  of  data 

Target  cell 
number 

SNR  [dB] 

0  5  10  15  20  25 

CR  theorem 

1.00  1.00  1.00  1.00  1.00  1.00 

48+53 

0,13 

Least  squares  method 

1.00  1.00  1.00  1.00  1.00  1.00 

100 

300 
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Example  12:  In  this  case,  the  signal  is  of  the  fonn 


n  ( 


and  only  the  Least  squares  method  can  be  used.  The  results  are  given  in  Table  3.8.12. 


Table  3.8.12  Summary  of  the  results  for  Example  12.  fs  =  —  fmax  with  5  signal  frequencies. 


Methods 

Error  1 

(time  domain) 

Error2 

(frequency  domain) 

Number  of 
data 

Target  cell 
number 

Least  squares  method 

- 

0.0000 

too 

300 

Example  13:  The  signal  is  given  by 

f(xk)=  I0e2xi42xk  +  ge2rd'2'6Xk  +  ge2Ml6Sxt  +  +6e2”'42-'r‘  +  Noise 

and  the  noise  is  characterized  as  in  (3-8-3)  and  the  SOI  is  same  as  in  Example  3.  The 
results  from  100  trials  are  described  in  Table  3.8.13. 


Table  3.8.13  Summary  of  the  results  for  Example  13.  fs  =  —  /max  with  signal  frequencies  and 

AWGN. 


Methods 

Probability  of  detection 

Number  of  data 

Target  cell 
number 

SNR  [dB] 

05  10  15  20  25 

Least  squares  method 

0.81  0.97  1.00  1.00  1.00  1.00 

100 

300 

Example  14:  The  signal  has  5  frequency  components  including  Gaussian  noise  and 
ground  clutter  has  been  added  so  that 


/(**)  =  i0g2 


+  8e 


+  9e 


+  2e 2 


+  6e 


+  noise  +  clutter 


The  noise  and  the  ground  clutter  are  the  same  as  described  by  (3-8-3)  and  (3-8-4)  and 
other  conditions  are  similar  to  the  previous  example.  The  simulation  result  is  given  in 
Table  3.8.14. 

In  summary,  when  the  sampling  frequency  is  much  less  than  the  highest  frequency  of 
the  signal,  only  the  Least  squares  method  and  the  CR  theorem  approach  can  be  applied  to 
compute  the  spectrum  of  a  nonuniformly  sampled  data  sequence.  In  the  presence  of 
ground  clutter  only  the  Least  squares  method  can  be  used.  The  performance  enhancement 
in  terms  of  the  maximum  resolvable  Doppler  frequency  is  given  by  (2-11,  2-14). 
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Table  3.8.14  Summary  of  the  results  for  Example  14.  fs  =  —  fmm 


with  5  signal  frequencies  with 


added  AWGN  and  ground  clutter  added. 


Methods 

Probability  of  detection 

Number  of 
data 

Target  cell 
number 

SNR  [dB] 

0  5  10  15  20  25 

Least  squares  method 

0.75  0.80  0.99  1.00  1.00  1.00 

100 

300 

3.8.4  Comparison  between  the  Least  Squares  Method  and  the  FFT 

From  the  results  of  the  previous  section,  when  the  Doppler  increment  due  to  a  multiple 
PRF  system  needs  to  be  large,  that  is  the  maximum  Doppler  exceeds  the  average  sampling 
rate,  only  the  clustering  algorithm  based  on  the  CR  theorem  and  the  Least  squares  method 
can  be  used.  The  results  from  both  of  the  methods  are  compared  in  tenns  of  stability  to 
noise  in  this  section.  The  quality  of  performance  depends  on  specific  situations.  We 
consider  a  test  signal  which  is  located  at  10kHz  with  unit  amplitude  sampled  by  4  PRFs  of 
1kHz,  1.35kHz,  1.6kHz  and  1.7kHz.  Gaussian  white  noise  having  real  and  imaginary 
components  with  equal  power  density  has  been  added  to  the  signal  as  in  equation  (3.8.3). 
The  signal  has  also  been  sequentially  sampled  at  32  sampling  points  by  using  all  the  PRFs. 
We  assume  that  an  error  occurred  if  the  amplitude  of  the  sidelobe  is  larger  in  magnitude 
than  that  of  the  actual  the  signal.  The  result  shows  that  the  Least  squares  approach 
outperformed  the  method  based  on  the  CR  theorem  and  the  FFT  method.  Figure  3.8.4  is  a 
plot  of  the  increase  in  the  noise  power  versus  the  probability  of  false  alarm  for  both  the 
methods  and  an  error  occurs  when  the  maximum  does  not  occur  at  the  real  signal 
frequency.  1000  simulations  have  been  carried  out  to  get  the  estimate  of  the  probability 
density  function.  Around  -2.5  dB  of  SNR,  the  CR  theorem  starts  to  fail  and  around  -5  dB 
of  SNR,  the  Least  squares  method  starts  to  fail.  Note  that  only  one  false  alarm  of  any  of 
the  4  PRF  may  lead  to  a  failure  of  the  method  based  on  the  CR  theorem  (clustering 
algorithm)  method.  Therefore  the  Least  squares  method  outperformed  the  method  based 
on  the  CR  theorem  by  2.5dB  of  SNR.  If  we  consider  a  5  or  6  PRF  system,  the  Least 
squares  method  will  provide  better  results  than  the  method  based  on  the  CR  theorem. 


3.8.5  Operation  Count 

The  speed  of  computation  associated  with  each  algorithm  is  also  an  important  factor. 
To  obtain  a  response  at  a  single  frequency  cell,  a  FFT  based  method  needs  logA 
multiplications  where  N  is  the  number  of  data  samples.  Table  3.8.15  shows  the  operation 
count  for  each  method.  A  filtering  operation  needs  2  FFT  and  one  vector  multiplication 
and  requires  2MogA  +  N  computation  for  the  evaluation  of  the  frequency  domain 
response.  A  matrix  inversion  needs  of  the  order  of  A3  operations  and  the  Hilbert 
transfonn  needs  2MogA  +  N  operations  to  evaluate  the  frequency  spectrum.  The  method 
based  on  the  CR  theorem  is  the  fastest  algorithm  and  the  Least  squares  method  is  also 
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within  an  acceptable  range.  The  other  techniques  require  0( N 3 )  operations  and  hence  the 
computational  complexity  becomes  large  when  N  becomes  a  large  number. 


Figure  3.8.4  Comparison  of  the  results  between  the  Least  squares  method  and  the  method  based  on  the  CR 
theorem  and  FFT.  A  false  alarm  occurs  when  the  sidelobe  level  is  greater  than  the  mainlobe.  As  SNR 
decreases  the  Least  squares  method  outperforms  the  FFT  method  by  about  2.5  dB. 


Table  3.8.15  Operation  count  of  each  method. 


Methods 

Operation  count 

From  equation 

Lagrange  polynomial 

0(N3) 

(3-1-1),  FFT 

Cauchy’s  method 

0(N3) 

Singular  value  decomposition. 

CR  theorem 

0(N  log  N) 

(3-2-28),  FFT 

Least  squares  method 

0(N2) 

(3-3-17),  Hilbert  transform 

QMF  method 

0(N3) 

Matrix  inversion,  filtering 

Iterative  method 

0(N3)x  Number  of  iterations 

(3-5-8),  (3-1-1),  filtering 

Orthogonal  polynomial  approach 

0(N3) 

(3-6-8),  Matrix  inversion 

Analog  frequency  approach 

0(N3) 

(3-7-6),  Matrix  inversion 

Chapter  4:  Conclusions 

This  report  addresses  the  problem  of  estimating  the  spectrum  from  a  set  of 
nonuniformly  spaced  data  for  applications  in  a  multiple  PRF  radar  system.  The  benefits  of 
using  a  direct  spectrum  analysis  to  a  set  of  nonuniformly  spaced  data  instead  of  using  a 
FFT  applied  to  the  method  based  on  the  CR  theorem  has  been  outlined.  The  direct 
spectrum  analysis  can  detect  multiple  targets  quite  easily  in  the  presence  of  the  ground 
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clutter  as  this  is  the  only  method  that  can  handle  signals  embedded  in  interferences  which 
are  of  similar  strengths  as  the  target. 

To  obtain  the  spectrum  from  nonunifonnly  sampled  data  various  methods  have  been 
studied  and  compared  with  each  other.  The  presented  methods  are: 

•  Polynomial  interpolation  (Lagrange  and  Cauchy  type) 

•  Chinese  remainder  (CR)  theorem  along  with  a  clustering  algorithm 

•  Least  squares  curve  fitting  of  a  nonunifonnly  spaced  complex  sequence 

•  Multi-resolution  analysis  (QMF  analysis) 

•  Iterative  method 

•  Orthogonal  polynomial  expansion  (Legendre  and  Hennite  polynomials) 

•  Estimation  of  the  analog  frequency 

The  numerical  simulations  have  been  performed  in  terms  of  different  frequency 
components  of  the  signal  for  each  method.  When  the  frequency  of  the  signal  is  smaller 
than  the  sampling  frequency,  most  of  the  methods  presented  including  the  interpolation 
methods  yield  acceptable  results.  As  the  frequencies  in  the  signal  come  close  to  the 
sampling  frequency,  most  of  the  time  domain  interpolating  techniques  fails.  The  Least 
squares  method,  the  method  based  on  the  CR  theorem  and  the  analog  frequency  technique 
can  successfully  analyze  signals  which  are  close  to  the  sampling  frequency.  If  the  signal 
frequency  is  higher  than  the  average  sampling  rate,  only  the  Least  squares  method  and  the 
method  based  on  the  CR  theorem  approach  can  be  used  for  the  analysis  of  the  data.  If 
there  exists  ground  clutter  along  with  multiple  signals,  only  the  Least  squares  method  can 
be  used  to  analyze  the  nonunifonnly  spaced  signal. 

Except  for  the  time  domain  interpolation  technique,  the  polynomial  interpolation 
method  and  the  orthogonal  polynomial  expansion  method,  all  other  methods  are  quite 
robust  to  AWGN  type  of  noise. 

Table  4.1  summarizes  the  characteristics  of  all  the  methods.  Here,  the  number  of  time 
domain  data  equals  N.  The  number  of  frequency  domain  samples  is  also  equal  to  N.  The 
number  of  PRF  equals  P  and  all  the  PRFs  are  considered  to  be  relatively  prime  numbers 
for  simplicity.  Performance  enhancement  is  given  in  terms  of  equation  (1-4)  when 
compared  to  a  single  PRF  system.  The  third  and  forth  columns  in  Table  4.1  represent 
whether  a  given  method  can  detect  a  target  successfully  in  the  presence  of  interferences 
and/or  ground  clutter. 

The  following  areas  of  research  may  be  of  interest  in  pursuing  in  the  future.  They  are 
listed  as  follows: 

•  Combination  of  the  various  algorithms  presented  in  this  report  for  enhancing  the 
estimation  of  the  Doppler  frequency. 

•  Extension  of  the  evaluation  of  the  spectrum  in  the  two-dimensional  Doppler  and 
Range  domain  by  using  nonunifonnly  sampled  data. 

•  Use  of  an  adaptive  filter  algorithm 

•  Simulation  using  real  data  from  a  practical  radar  system  which  often  encounters 
non-  stationary  or  non-homogeneous  environments. 
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Table  4. 1  Summary  of  the  characteristic  of  each  method 


Method 

Operation 

count 

Performance 

enhancement 

Multiple 

interferences 

Burst  in  the 
clutter  band 

Polynomial  interpolation 
method 

0(N3) 

,  I  PRFi 

8PRFref 

Yes 

Yes 

Least  square  method 

0(N2) 

l.c.m.  (PRFj ) 

PRFref 

Yes 

Yes 

QMF  method 

0(N 3) 

UPRFi 

PRFref 

Yes 

Yes 

Orthogonal  Expansion 

0(N3) 

.  I  PRFt 

8  PRFref 

Yes 

Yes 

Analog  frequency  method 

0(N3) 

YjPRFj 

PRFref 

Yes 

Yes 

CR  theorem 

0(N  log  N) 

l.c.m.  ( PRFi ) 

PRFref 

No 

No 

Bibliography 


[1]  M.  I.  Skolnik,  Radar  Handbook,  1970,  McGraw-Hill  Book  Company,  New  York. 

[2]  H.  R.  Raemer,  Radar  System  Principles,  1996,  CRC  press,  Boca  Raton,  FL. 

[3]  A.  Ludloff  and  M  Minker,  “Reliability  of  Velocity  Measurement  by  MTD  Radar,” 
IEEE  Transactions  on  Aerospace  and  Electronic  Systems,  AES-23,  4,  July  1985,  pp.  522- 
528. 

[4]  I.  Vrana,  “Optimum  Statistical  Estimates  in  Conditions  of  Ambiguity,”  IEEE 
Transactions  on  Information  Theory,  Vol.  39,  No.  3,  May  1993,  pp.  1023-1030. 

[5]  G.  Trunk  and  S  Brockett,  “Range  and  Velocity  Ambiguity  Resolution,”  IEEE  National 
Radar  Conference,  1993,  pp.  146-149. 

[6]  O.  Ore,  Number  Theory  and  Its  History,  1949,  McGraw-Hill  Book  Company,  New 
York. 

[7]  Y.  Hua  and  T.  K.  Sarkar,  “Matrix  pencil  and  system  poles,”  Signal  Processing,  Vol. 

21,  No.  2,  October  1990,  pp.  195-198. 

[8]  T.  K.  Sarkar.  and  O.  Pereira,  “Using  the  Matrix  Pencil  Method  to  Estimate  the 
Parameters  of  a  Sum  of  Complex  Exponentials,”  IEEE  Antenna  and  Propagation 
Magazine,  Vol.  37-1,  pp.  48-55,  Feb  1995. 

[9]  P.  S.  Naidu,  Modern  spectrum  analysis  of  Time  Series,  1996,  CRC  press  Boca  Raton, 
FL. 

[10]  P.  Vanicek,  “Further  Development  and  Properties  of  the  Spectral  Analysis  by  Least 
squares”,  Ap.  Space  Sci.,  Vol.  4,  pp.  10-33,  1971. 


72 


[11]  N.  R.  Lomb,  “Least  squares  frequency  analysis  unequally  spaced  data”,  Airophysics 
and  space  science,  Vol.  39,  pp.  447-462.  May.  1975. 

[12]  F.  J.  Beutler,  Error-Free  Recovery  of  Signals  from  Irregularly  Spaced  Samples, 
SIAM  Review,  Vol.  8,  No.  3,  pp.  328-335,  July  1966. 

[13]  H.  S.  Black,  Modulation  Theory,  1953,  Van  Nostrand,  New  York. 

[14]  M.  Abramowitz  and  I.  A.  Stegun,  Handbook  of  mathematical  functions,  1970,  Dover, 
New  York. 

[15]  I.  Bilinskis  and  A.  Mikelsons,  Randomized  Signal  Processing,  1992,  Prentice  Hall, 
New  York. 

[16]  R.  Adve  and  T.  K.  Sarkar,  The  effect  of  noise  in  the  data  on  the  Cauchy  method. 
Microwave  and  Optical  Technology  Letters,  Vol.  7,  No.  5,  April  5  1994,  pp.  242-247. 

[17]  K.  Kottapalli,  T.  K.  Sarkar,  Y.  Hua,  E.  Miller  and  G.  J.  Burke,  “Accurate 
Computation  of  Wide-band  Response  of  Electromagnetic  Systems  Utilizing  Narrow-band 
Infonnation”,  IEEE  Trans.  On  MTT,  Vol.  39,  pp.  682-688,  April  1991. 

[18]  S.  V.  Huffel,  “Analysis  of  the  Total  Least  Squares  Problem  and  its  use  in  Parameter 
Estimation”,  Ph.  D.  Dissertation,  Department  Electro technick,  Kattrolicke  Universiterit 
Leuven. 

[19]  J.  D.  Scargle,  “Studies  in  astronomical  time  series  analysis.  II.  Statistical  Aspects  of 
spectral  analysis  of  unevenly  spaced  data”,  The  Airophysical  Journal,  Vol.  263,  pp.  835- 
853,  Dec.  1982. 

[20]  J.  H.  Home,  and  S.  L.  Baliunas,  “A  prescription  for  period  analysis  of  unevenly 
sampled  time  series”,  The  astrophysical  journal,  Vol.  302,  pp.  757-763,  Mar.  1986. 

[21]  S.  F.  Mello,  “Estimation  of  periods  from  unequally  spaced  observations”,  Astron.  J, 
Vol.  86,  No.  4,  pp.  619-624,  Apr.  1981. 

[22]  D.  D.  Meisel,  “Fourier  transforms  of  data  sampled  at  unequal  observational 
intervals”,  Astro.  J,  Vol.  83,  No.  5,  pp.  538-545,  May.  1978. 

[23]  W.  H.  Press,  S.  A.  Teukolsky,  W.  T.  Vetterling,  and  B.  P.  Flannery,  Numerical 
recipes  in  c.  Cambridge:  University  press,  1992. 

[24]  P.  A.  Gorry,  “General  least  squares  smoothing  and  differentiation  of  nonuniformly 
spaced  data  by  the  convolution  method”,  Anal.  Chem.,  Vol.  63,  pp.  534-536,  1991. 

[25]  P.  P.  Vaidyanathan,  Multirate  Systems  and  Filter  Banks,  1993,  Prentice-Hall, 
Englewood  Cliffs,  NJ. 

[26]  P.  P.  Vaidyanathan,  and  V.  C.  Liu,  “Classical  sampling  theorems  in  the  context  of 
multirate  and  polyphase  digital  filter  bank  structure”,  IEEE  Trans,  on  Acoustics,  Speech 
and  Signal  Proc.,  Vol.  ASSP-36,  pp.  1480-1495,  Sept  1988. 

[27]  P.  P.  Vaidyanathan,  and  V.  C.  Liu.  "Effective  reconstruction  of  band-limited 
sequences  from  nonuniformly  decimated  versions  by  use  of  polyphase  filter  banks”,  IEEE 
Trans,  on  Acoustics,  Speech  and  Signal  Proc.,  Vol.  38,  pp.  1927-1936,  Nov  1990. 

[28]  I.  W.  Sandberg,  “On  the  properties  of  some  systems  that  distort  signals-I”,  Bell  Syst. 
Tech.  J.,  Vol.  42  ,  pp.  2033-2047,  Sept.  1963. 

[29]  R.  G.  Willy,  “Recovery  of  band-limited  signals  from  unevenly  spaced  samples”, 
IEEE  Trans.  Cominun.,  Vol.  COM-26,  No.  1,  Jan  1978. 

[30]  F.  Marvesti,  M.  Analoui,  and  M.  Gamshadzai,  “Recovery  of  signals  from  nonuniform 
samples  using  iterative  methods”,  IEEE  Trans.  Signal  Processing,  Vol.  39,  pp.  872-878, 
Apr.  1991. 


73 


[31]  I.  W.  Sandberg,  “The  reconstruction  of  band-limited  signals  from  nonuniformly 
spaced  samples”,  IEEE  Trans.  Circuits  and  Systems,  Vol.  41,  pp.  64-66,  Jan.  1994. 

[32]  Y.  Park,  and  M.  Soumekh,  “Reconstruction  from  unevenly  spaced  sampled  data  using 
iterative  methods”,  IEEE  Trans.  Signal  Processing,  Vol.  43,  pp.  303-308,  Jan.  1995. 

[33]  M.  M.  Rao,  T.  K.  Sarkar  and  R.  S.  Adve,  “Simultaneous  extrapolation  in  time  and 
frequency  domains  using  Hermite  expansions”,  IEEE  Trans.  Antenna  sand  Propagat.,  Vol. 

[34]  Y.  C.  Jenq,  “Perfect  Reconstruction  of  Digital  Spectrum  from  Nonuniformly  Sampled 
Signals”,  IEEE  Trans.  Instrument  and  Measurement,  Vol.  46,  No.  3,  pp  649-652,  June 
1997. 

[35]  Y.  C.  Jenq,  “Digital  Spectra  of  Nonuniformly  Sampled  signals:  Fundamental  and 
High  Speed  Waveform  Digitizers”,  IEEE  Trans.  Instrument  and  Measurement,  Vol.  37, 
pp.  245-251,  June  1988. 

[36]  Y.  C.  Jenq,  “Digital  Spectra  of  Nonuniformly  Sampled  signals:  Digital  look  up 
Tunable  Sinusoidal  Oscillators”,  IEEE  Trans.  Instrument  and  Measurement,  Vol.  37,  pp. 
358-362,  Sept  1988. 


74 


Appendix  A:  The  Sampling  theorem  for  a  randomly  sampled 
data 

If  the  sampling  of  the  data  were  to  be  completely  random,  and  a  Poisson  process,  then 
the  spectrum  would  be  completely  free  from  aliasing.  This  can  be  shown  from  the 
development  of  the  Fourier  transform  as  applied  to  a  uniformly  spaced  discrete  sequence. 
Consider  a  uniformly  sampled  signal  x(nT )  as 


oo 

x{nT )  =  x(t)  ^S(t  -nT), 

n=— oo 


(A-l) 


Then  its  Fourier  transform  is  given  by 


oo 

X(oj)=  Jx(«r)  e-jmTdt  =  Zx("r) e 


—jconT 


(A-2) 


The  Fourier  transfonn  is  periodic  since 


X\ 


In 
oo  +  — 

V  T  J 


i  2;r  |  _ 

-/  co-\ - nT 

1  T 


=  £x(nT)  e-^T-J2m  =  ei<mT  =X{co). 

n=— oo  n=— oo 

(A-3) 

For  the  case  of  the  nonuniformly  sampled  data,  the  Fourier  transform  will  be  defined  from 


00  oo 

X{po)=  Jx(0  e~J(0,"dt  =  ^x(c)  e~JC<  .  (A-4) 

-oo  «=-oo 

Here,  X(ru)  is  periodic  only  when  the  exponential  part  of  X(co)  is  periodic  as  shown  in 
the  uniformly  sampled  case.  That  is,  eJ03t'"  =  should  hold  for  all  tn .  If  the  time 

steps  are  not  detenninistic  quantities  then  the  probability  of  tm  and  tnj  having  the  same 
2  71 

time  increments,  — ,  as  for  the  uniformly  sampled  case,  for  the  entire  time  duration  of  the 
co 

signal  will  be  zero.  Therefore  X(<y)  is  not  periodic  in  co  and  x(t)  does  not  need  to  be  a 
band-limited  sequence  that  leads  to  an  alias  free  condition. 
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Appendix  B:  Matrix  Pencil  Method  (MPM) 


MPM  is  a  method  to  fit  a  uniformly  spaced  data  sequence  by  a  sum  of  complex 
exponential.  Sarkar  and  Hua  [7,  8]  described  in  details  this  method.  It  is  summarized  here 

for  completeness.  The  sampled  signal  x[kTs )  is  to  be  modeled  by  a  sum  of  complex 
exponentials,  i.e., 

M  M 

x[kT)~YR:c'kl  (B-l) 

l'=l  i=\ 

where  R:  =  Residues  or  complex  amplitudes, 

=-at  +  jcol , 
at  =  Damping  factors, 
col  =  Angular  frequencies, 
es‘T*  =  z.  for  i  =  1 ,  2,  . . . ,  M. 

The  objective  is  to  find  the  best  estimates  of  M,  Ri  and  z;  from  x(kTsj . 

We  can  define  matrices  [ij]  and  [k,]  as  follows  (Assume  we  have  N  sampled  data 
points): 


x{0) 

x(l) 

•  x(L  -  1)  “ 

M= 

x(l) 

x(2) 

x(L) 

A 

l  . 

l 

1)  x(N-L)  •• 

M 

i= 1 

M 

Zra 

i=\ 

M 

■  ZR^L~ 

i= 1 

= 

M 

Z  Rizi 

i=\ 

M 

z  =  l 

M 

-  Zral 

i=l 

M 

_  i=l 

1 

Jnj 

i 

M 

-  Zra”~ 

i= 1 

x(l) 

x(2) 

x(Z) 

[*,]= 

x(2) 

x(3) 

•  x(L  + 1) 

x(N-L) 

x(A-Z  +  l)  •• 

■  4^-1)! 
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M 

.■—1 

M 

i>a2  ■ 

M 

■  I >,V 

; _ 1 

i=i 

M 

Zra2 

i= 1 

i=\ 

M 

i>,  v  • 

i= 1 

1  =  1 

M 

i= 1 

M 

D  N-L 

LR:zi 

i= 1 

M 

V  n  N-L+l 

LRizi 

i= 1 

JN 

1 _ 

(B-3) 


where  L  is  called  the  pencil  parameter.  L  is  chosen  in  between  N/3  to  N/2  for  efficient 
noise  filtering  [7]. 

Since  we  do  not  know  how  many  frequency  components  exist  in  the  signal,  the  number 
of  estimated  frequencies  M  should  be  determined  using  some  criteria.  Typically  the 
singular  values  beyond  M  are  set  equal  to  zero.  The  way  M  is  chosen  is  as  follows  [8], 
Consider  the  singular  value  ac  such  that 

(B-4) 


where  p  is  the  number  of  significant  decimal  digits  in  the  data.  One  can  write 


LA]  = 


fJV-I-1)  JN-L-I) 


(N-L-l) 


Consider  the  matrix  pencil 


[K,]  =  [z,JifJ[Z(1][z,], 

1  1-1 

Z1  Z2  ' "  ZM 


1  Zj  •••  z[L~x) 

1  z2  •••  z*i_1) 

[z.1- 

}  Z  M  •••  ZM_1)_ 

c/h/y[z|  ,z2 , •  •  • ,  zM  | , 

R0]=diag[Rl,R2,---,RM]. 

(B-5) 

(B-6) 

(B-7) 

( N-L)xM 

(B-8) 

(B-9) 

(B-10) 
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(B-l  1) 


[^]-A^\=\z>\Ri\z«\-x wim- 

Therefore  A  =  z;. ,  for  i  =  1,  2,  M  would  be  the  exponentials  determined  from  the 
generalized  eigenvalue  problem. 


[rM  W- 


(B-12) 


It  can  be  shown  that  this  is  the  same  as  solving  the  ordinary  eigenvalue  problem 

frrfc]}-*  w.  (B-13) 


where  Y]  1  is  the  Moore-Penrose  pseudoinverse  of  Yl  which  is  defined  by 


[irMfomn*;]'. 


(B-l  4) 


Once  A,  =  z.  are  known,  the  residues  Rt  are  solved  for  from  the  following  least  square 
problem 


x(0) 

1 

1 

1 

"*ll 

x(l) 

= 

Z1  Z2  ' 

••  ZM 

r2 

(B-15) 

_x{N  -  l)j 

ZN-X 

L  i 

ZN -1  . 
z2 

...  ZN -1 
^  M  J 

Rm1 

The  frequency  component  is  computed  from 

(Oj  =  Im[ln(z; )]  (B-l  6) 

and  the  magnitude  A  t  for  a  single  frequency  co ;  is  evaluated  from 

Aiej(0i  =  Rie~a‘  +  Jm‘  =  Ri[-Re{zi)\  (B-17) 
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Appendix  C:  Proof  of  (3-5-6) 

Let  A  =  PQ  and  note  that 

Re(  Af  -Ag,f-g)  =  Re(gf  -  Qg,  Pf  -  Pg)  =  Refe/'  -Qg,f-g)  =  >k,  || /  -  gf  (C-l) 

for  all  /,ge  k  since  P  is  self-adjoint  operator. 

The  equation  P[  =  Af  is  equivalent  to  /  =  Af  ,  where  A  f  =  ch  +  f  -cAf  and  c  is 
any  nonzero  constant.  The  following  calculation  shows  that  A  is  a  mapping  from  k  into 
k.  It  is  a  contraction  when  c  =  kt(k2^j  : 

Af-Ag  =  \\f  -  g  -  cAf  +  cAgf  =\f  -  g\ -2c^e(Af  -  Ag,f  -  g)  +  c2\Af  -  Ag\ 


<  (l  -  2ckx  +  c2k2  )||/  -  g||2 ,  c>  0.  (C-2) 

So  k\2  <  k2 ,  when  (1  -  2cki  +  c2k2 )  >0  for  all  c>  0 . 

2  f  (  k2^ 

Af-Ag  <  l--H||/-g||\o<  (C-3) 

V  k 2  J  V  & 2  ' 

The  Schwartz  inequality  shows  that  the  last  equation  can  be  stated  as 

I  Ik  -  Ag\-\\f  -  g||  >  |  (Af  -  AgJ  -  g)|  >  k,  1/  -gf.  (C-4) 

Thus 

\\Af-Ag\>k\f-g\.  (C-5) 

With  /  =  A~\  and  g  =  A ~lh2 , 

| h{  -h2\\>  k^A~\  -A~lh2f .  (C-6) 
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Appendix  D:  Proof  of  Uniqueness  of  the  Solution 

Let  f,g  etc  and  let  Q  be  a  mapping  of  k  into  H  such  that  ( Qf  -  Qg,f  -  g)  vanishes 
only  if  /  =  g .  Then  if  the  equation  h  =  PQz  has  a  solution  z  e  /r ,  it  is  unique. 

Proof:  Assume  that  PQzx  =  PQz2  where  z,  ,z2  e  k  .  Since  P  is  a  self-adjoint  operator, 

(Qz i  -Qz2-zx  -z2)  =  (Qz1  -Qz2,Pz1  -  Pz2)  =  (pQzx  -PQz2,z1  -z2)  =  0.  (D-l) 

Hence  zt  =  z2 . 

Appendix  E:  Proof  of  (3-5-12) 

The  left-hand  side  of  (3-5-12)  can  be  written  as 

p[xk(t)  -  xk_x(t)\-\P<^xk(t)  -  xk_x(t)\ |  =|xt(f)-xt_j(f)|r  +L2  P<^xk{t)-xk_x(t)\ | 

-  2Lj  P[xk{t)  -  xk_x{t)]PQ[xk{t)  -  xk_x[t)]dt .  (E-l) 

We  want  to  show  that  there  exist  positive  real  numbers  kx  and  k2  such  that 

J  p(xk  (0  ~  xk- 1  (t))pQ{xk  0)  -  xk_t  (t t))dt  >  kx  ||x,  (0  -  xk_x  (Of  (E-2) 

and 

\pQ{xk{t)-xk_x{t))  <A:2||x,(0-W-i(0|2.  (E-3) 

To  prove  (E-2),  we  first  note  that 

J  p(xk  (0  -  *k- 1  (0 )pQ{xk  (0  -  W-i  (t))dt  =  J  (x,  (0  -  xk_x  (t))Q(xk  (0  -  xk_x  (i t))dt ,  (E-4) 

because 

P 2  =  P, 

Pxk  =  xk  ■ 

Since  the  operator  P  is  self-adjoint  (i.e.,  J  ( Px)ydt  =  J  x(  Py)clt ),  (E-4)  can  be  rewritten  as 
J  p{xk{t)-xk_x{t))pQ{xk(t)-xk_x(t))dt  =  S[w(h)-w-i(h  )]"•  (E-5) 
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For  a  band  limited  signal  x(t)  one  can  find  positive  numbers  A  and  B  which  satisfy 


A<—, - ^<B 


Therefore 


J  p(xk  (*)  -  xk_x  ( t))PQ(xk  {t)  -  xk_x  ( t))dt  >  A\xk  (t)  -  xk_x  (f  )f 


and  kx  =A  exist. 

To  prove  (E-3),  one  can  write 


|  po(xk  (0  -  xu- 1  (0)f  =  J  pQ(xk  0)  -  xk~  i  (0)Z  ( xk  (0  -  xk- 1  (*))§(*  -  h  )dt 

i 

=  Z  J  pQ(xk  (0  -  xk- 1  W)(w  )  -  xk- 1  (h  ))5(^  -  h  )dt 

i 

=  H[pQ(.xk  (0  -  **-i  (0)]*  =  h  (**  (h )  -  **-i  (h )) 

2iXr 

Z  k  k )  -  x*-i  k  )|  Z  \pQ(xk  (0 

.  i  J  L  i 

<  B\\xk  {t)~  xk_x ( t)\  ■  PQ(xk (t)  -  xk_x ( t)) 


k-WJ  I  t=u 


(E-6) 


(E-7) 


(E-8) 


We  have 


\pQ{xk(t)  -  xk_x{tj)  <  . 

(E-9) 

Therefore  k2  - 

=B  exist.  From  (E-l)-(E-3),  we  get 

p[xk{t)-xk_x{t)\-XP^xk{t)~xk_x{t)^ 

(E-10) 

and 

r  =  1  +  A  2k2  -  2  A  kx . 

(E- 11) 

In  order  to  satisfy  (3-5-12),  r  has  to  be  in  the  region  0  <  r  <  1 .  Given  a  particular  kx  and 
k2 ,  X  has  to  fall  within  the  following  region  of  convergence  for  the  iterative  relationship 
given  in  (3-5-8)  to  hold. 

2k, 

0<A<— (E-12) 
k2 


k{  <  k 


2  * 


(E-13) 
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LIST  OF  PRINCIPLE  SYMBOLS 


c  speed  of  light,  c  =  3  •  1 0 8  m/sec 

At  round  trip  transmit  time  of  the  wave  transmitted  and  reflected  back  to  the 
origin 

/'  frequency  of  the  transmitted  signal 

rmax  maximum  range  without  any  ambiguity 
Fmax  maximum  Doppler  velocity  without  any  ambiguity 

/max  maximum  Doppler  frequency 

Aq  wavelength  corresponding  to  the  carrier  frequency  of  the  radar 

a(t )  envelope  of  the  signal 

Fc  carrier  frequency 

'P(t)  phase  function 

v  target  moving  at  speed 

a  scale  factor  controlled  by  the  Doppler  effects 

Fd  Doppler  velocity  of  the  target 

T  width  of  the  pulse  in  a  period 

T  period  of  the  base-band  pulse 

A 

A {t,a)  complex  ambiguity  function  of  g(t) 

g(t,a )  received  signal 

PRFi  7-th  pulse  repetition  frequency 

ID  maximum  Doppler  increment 

IR  maximum  range  increment 

PRF2  pulse  repetition  frequency  for  comparison  to  single  PRF  system 

l.c.m  Least  common  multiplier 

g.c.d  greatest  common  divisor 

xk  sampling  points 

f[xk )  value  of  the  signal  at  xk 

f0  target  Doppler 

x0  target  range 

A  mat ,  R  mm  matrices  used  in  the  Cauchy’s  method 

p,  q  degrees  of  numerator  and  denominator  used  in  the  Cauchy’s  method 
c(j)  error  in  the  clustering  algorithm 

r  delay  parameter  enables  to  select  any  arbitrary  origin  of  time 

F  mean  square  difference  in  the  Least  square  method 

E  (co  )  spectrum  estimate  using  a  Least  square  method 
A  unknown  magnitudes  vector 

B  matrix  used  to  get  the  magnitude  vector  A 

U,  V  intennediate  steps  in  the  QMF  approach 
H  band-pass  filter 
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M 


total  number  of  time  steps  in  one  period  where  decimation  repeatedly 
occurs 

N  number  of  known  samples  in  the  M  time  intervals 
Q  band-limiting  operator  (low  pass  filter), 

P  ideal  nonuniform  sampling  operator 

X  coefficient  of  iterative  method 

Tk  uniform  sampling  point  with  average  interval 
ATk  k- th  interval  of  uniform  average  spacing 

hn  associate  Hennite  function  of  degree  n 
Hn  Hennite  function  of  degree  n 
/, ,  l2  scaling  factors 
Pn  Legendre  polynomial  of  degree  n 
rm  m- th  deviation  from  uniform  sampling 
Xa  (co)  analog  spectrum 
Xd  (<z>)  digital  spectrum 
Y  matrix  used  to  obtain  Xa  (co) 
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